Neural Network-Based Closed-Loop Deep Brain Stimulation for Modulation of Pathological Oscillation in Parkinson’s Disease

Aiming at the problem that the Proportional-Integral-Derivative (PID) control strategy needs to readjust controller parameters for different Parkinson’s disease (PD) states. This work proposes an improved control strategy that considers an artificial neural network control scheme. A backpropagation neural network (BPNN) controller is designed to solve the above problem and further to improve the performance of the closed-loop control strategy. The training data set of the BPNN controller is obtained by controlling eight different PD states (PDa - PDh) by the PID controller and the BPNN controller is trained by the training data set to obtain a set of optimal weights. By modulating other different PD states (e.g. PD1 - PD3), the effectiveness of the PID-structure controller and BPNN controller are compared. We find that the BPNN controller can modulate different PD states without changing the controller parameters and reduce energy expenditure by 58.26%. This work is helpful for the design of more effective closed-loop deep brain stimulation (DBS) systems for clinical applications and provides a framework for the further development of closed-loop DBS.


I. INTRODUCTION
Parkinson's disease (PD) is a neurodegenerative disease and its symptoms can be complicated to manage. The main motor symptoms of PD are loss of autonomic movement, decreased autonomous movement, muscle stiffness, and resting tremor of the limb [1]. The currently recognized origin of PD is the degradation of the substantia nigra pars compacta (SNc) dopaminergic neurons. However, the pathogenesis of PD is still unclear. It is generally believed that the Parkinsonian state may be related to abnormal synchronous oscillatory activity occurring in the basal ganglia-thalamocortical circuit [2]. The Parkinsonian state is associated with changes in neural activity, including altered neuron firing The associate editor coordinating the review of this manuscript and approving it for publication was Shuping He . rates, increased burst firings, and enhanced beta oscillatory activity across the basal ganglia-thalamo-cortical network [3], [4]. Therefore, these changes in neural firing activities are of great significance for guiding the research of PD. Deep brain stimulation (DBS) has the advantages of reversibility and adjustability, and is now increasingly favored by patients with refractory PD [5]- [7]. DBS, as a treatment approved by the US Food and Drug Administration [8], stimulates basal ganglia nuclei (such as Subthalamo nucleus (STN), Globus Pallidus interna (GPi) [9]) and the ventral intermediate nucleus (Vim) [10] in open-loop form using continuous high-frequency (greater than 100Hz) electrical stimulation to achieve modulation of the pathological neuronal rhythm of networks in patients with PD [11]. Clinical studies have shown that DBS can achieve a remarkable effect and has become the first choice for patients VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ with refractory PD. However, current clinical application mainly uses open-loop DBS, which has disadvantages. The parameters of DBS are set according to clinical experience, and the fixed stimulation pulse sequence can be difficult to adapt to differences between individual patients, requiring many visits. Additionally, continuous application of constant high-frequency electrical stimulation can result in excessive stimulation that increases the risk of side-effects and the need for surgery to replace the battery due to excessive energy consumption [12], [13].
Closed-loop DBS, as a treatment for PD, is an active area of research. Closed-loop DBS can operate on-demand and adjust the stimulation parameters in real time according to changes of the patient's physiological signals [14]- [16], which can effectively overcome the shortcomings of current open-loop DBS methods. The closed-loop DBS system is able to adjust the stimulation parameters in real time based on changes in the patients' clinical states and to apply specific stimuli parameters according to the needs of different patients [17], [18]. The DBS waveform can be realized based on a combination of control theory and real-time variations of the neural states. The effectiveness of closed-loop DBS has been successful demonstrated using computational models. Pyragas et al. previously proposed a PID control strategy [19]. They presented a demand-controlled method for desynchronization of globally coupled oscillatory networks utilizing a configuration with an observed and stimulated subsystem. The stimulated subsystem is subjected to a PID feedback derived from the mean field of the observed subsystem. Dunn et al. proposed PID closed-loop DBS control scheme [20]. However, a weakness in the PID control strategy is the need for readjusting the controller parameters for different PD states. That is, a group of PID parameters that work well for one Parkinsonian state do not work well for another Parkinsonian state. This requires re-adjusting the parameters of the PID controller to work effectively in different Parkinsonian states. In order to solve this problem, a backpropagation neural network (BPNN) controller is introduced in this work. The BPNN, also known as a feedforward neural network, has been well studied [21], [22]. Gharehchopogh et al. proposed two types of artificial neural networks to classify effective diagnosis of PD [23]. Ye et al. used artifcial neural networks to determine weights to scale the controller parameters. Their results demonstrated that Zishenpingchan granules reduced the occurrence of motor complications, and were useful for mitigating dyskinesia and non-motor symptoms of PD [24]. Considering the nonlinear mapping, self-learning and self-adaptive abilities and fault tolerance of BPNN, the present work introduces a BPNN controller to realize the control of DBS to the different PD states.
The main remainder of the paper is organized as follows: Section II describes a computational model of PD based on physiological structures within the basal ganglia -thalamocortical neural network. The BPNN closed-loop control scheme and evaluation criterion are also given in Section II. In Section III, simulation results are presented and the performance of the control is analyzed. For three different Parkinsonian states (PD1 -PD3), the BPNN-based closedloop DBS stimulation was applied to STN. Simulation results show that the BPNN controller achieves better control performance without changing the controller parameters, as compared with the PID algorithm. Finally, conclusions are given in Section IV.

A. MODELS OF BASAL GANGLIA-THALAMO-CORTICAL CIRCUIT
A schematic diagram of the basal ganglia-thalamo-cortical computational model is shown in Fig. 1. The network model constructed in this work includes six main types of neuron groups: pyramidal (PY), interneuron (IN), Globus Pallidus externa (GPe), GPi, STN and thalamus (TH). Each nucleus contains 100 neurons. The red arrows indicate excitatory synaptic connections and the blue arrows indicate inhibitory synaptic connections. The structure of the basal gangliathalamo-cortical circuit model adopts the sparse connections as in [25], which originated from Rubin and Terman's work [26].
However, to approach the biologically realistic, we make some modifications from Lu's model [25] in this work. The first modification is changing the model structure to a random pattern [27]. In our modified basal ganglia-thalamo-cortical computational model, each GPi neuron receives excitatory synaptic inputs from three randomly selected STN neurons and each TH neuron receives inhibitory synaptic inputs from one randomly selected GPi neuron in this model. In the modified cortical model, each PY neuron receives excitatory synaptic inputs from five randomly selected TH neurons and each IN neuron receives excitatory synaptic inputs from one randomly selected PY neuron.
All the neurons are simulated by the Izhikevich model, whose mathematical description is as follows [25]: with the reset equation after firing: where the variable V represents the membrane voltage and r represents the membrane recovery variable. In order to characterize the dynamics of the neurons in the basalthalamo-cortical neural network and effects of the DBS [26], the values of a, b, c, d and I are modified from [25] as shown in Tab. 1, which is the second modification of our model. The third modification is that our model consider the synaptic delay effects [27]. By introducing the synaptic delay, the synaptic current from the jth neuron to the ith neuron is changed to where τ is the delay time with values given in Tab. 1.
Besides, g ij describes the synaptic connection strength, and the changes of values of g ij reflect different states in Tab. 1 [28]. In (4), the synapse variable S j , as a fraction of the post-synaptically bound neurotransmitter, obeys first-order kinetics [25] and E syn is the reversal potential with E syn = 0mV for excitatory synapses and E syn = −80mV for inhibitory synapses. Based on these modifications of the structure and parameter, the DBS effects can be reflected in the neuron firing activity variations of all the nuclei especially the cortex. Thus, the modified basal ganglia-thalamo-cortical closed-loop neural network is more appropriate to be used in the exploration of closed-loop DBS optimization.

B. DESIGN OF THE CLOSED-LOOP DBS CONTROL STRATEGIES
The closed-loop DBS control system is constructed to modulate the behavior of the basal-thalamo-cortical neural network. A block diagram of the closed-loop control system constructed in this work is shown in Fig. 2. The switch K 1 and K 2 are used to choose the controller. When K 1 is switched on, PID controller is applied, and when K 2 is switched on, BPNN controller is applied.
Cortical information has great importance to reflect the DBS effects for PD patients and is much easier to be recorded compared with the neural signals in the deep brain. For better integration with EEG signals in the future, the average membrane potential of the PY nuclei neurons is selected as the feedback signal. In the closed-loop control system block diagram shown in Fig. 2, the controller is PID controller or BPNN controller, DBS is equivalent to an actuator, basal ganglia-thalamo-cortical network acts as the controlled object and y PY is the controlled variable. y PY is the average membrane potential signal of the PY neurons in the closed-loop DBS control system. y PY is the average membrane potential of the PY neurons in the normal state and used as the expected signal, u is the controller output signal, and z is the stimulus signal added to the controlled object after modulation by DBS, and e is the error between y PY and y PY . The average membrane potential for nucleus I can be caluculated by (5).
An example of a DBS waveform is included in Fig. 2. The magnitude of the DBS is m, and in closed-loop pattern, the stimulus amplitude m is equal to control signal u. The stimulated nucleus is the STN based on a typical clinical application. A phenomenon observed in primate experiments is that a DBS stimulus with a frequency higher than 100Hz acting on the STN can relieve PD symptoms, and is usually ineffective when the frequency is lower than 50 Hz [29], [30].
In an experimental rat model, the effective frequency range is 50-130 Hz, and the therapeutic effect is saturated when the frequency is greater than 130 Hz [31]. The standard DBS stimulation signal consists of a sequence of high-frequency VOLUME 8, 2020 biphasic electrical pulses above 100 Hz, typically at 130 Hz [27].
Thus, in this work, the frequency of the pulse sequence is chosen as 130Hz. In order to avoid tissue damage, a set of phase-inverted electrical pulses with longer duration and lower amplitude are provided to balance the net injected charge [32]. Thus, biphasic charge-balanced DBS is considered as shown in Fig. 2. The short first pulse with duration 300µs (between 30µs to 450µs as suggested in [33]) is followed by a charge-balancing second phase of opposite polarity with ten times longer duration than that of the short first phase pulse [34]. The amplitudes of the two phases are correspondingly adjusted, such that the net charge is zero [35].

1) PID ALGORITHM FOR THE CLOSED-LOOP CONTROLLER
When K 1 is switched on, PID controller is applied. Using a classical PID controller, the control signal u can be represented by u PID , which is adjusted as follows: where K p is the proportional coefficient, K i is the integral coefficient, and K d is the differential coefficient.  times individually. At the same time, record the difference between the average membrane potential of the PY nucleus in each PD state (y PYay PYh ) and the normal state (y PY ). Therefore, eight sets of corresponding input errors are obtained, as shown in e a (n)e h (n) in Fig. 3. By repeatedly adjusting the parameters of the PID controller, a corresponding set of better control signals are obtained in each PD state, and as shown in u a (n) − u h (n) in Fig. 3. Therefore, eight sets of corresponding input data (e a (n)e h (n)) and eight sets of corresponding output data (u a (n) − u h (n)) are obtained. The simulation time for each PD state is 5 seconds, and the simulation step is 0.1 milliseconds, so each PD state group includes 50,000 samples of data. Because the amount of data is too large, downsampling is used and one data sample is taken every 10 steps. There are a total of 8 different PD states so obtain an 8 × 5000 = 40000 sample input-output data set. The data was normalized and disordered before training the neural network, and then the BPNN was trained.

b: TRAINING ALGORITHM
This work adopts the error back propagation algorithm. The BPNN controller uses the error correction learning method to continuously modify and optimize the connection weights according to the expected output and the obtained actual output error to achieve the desired mapping relationship.
As shown in Fig. 2, the BPNN controller is divided into three layers: the input layer, the hidden layer, and the output layer. The number of neurons in the input layer and the output layer is determined by specific questions. In this work, the number of neurons in the input and output layers is 1, and the number of neurons in the hidden layer is 5 (discussed in section III). The input and output of the network are e and u, respectively, and w ij and v jk are the connection weights. The activation function of the hidden layer uses the sigmoid function. The flowchart of the BPNN training process is shown in Fig. 4. The training process of BPNN can be divided into the following two stages: a) The forward phase of the working signal. The calculation formula for the jth neuron input of the hidden layer is as shown in (7): l j (n) = w 1j (n)e(n) j = 1, 2, . . . , 5 The result of (8) is the output of the j th neuron of the hidden layer: h j (n) = Sigmoid l j (n) j = 1, 2, . . . , 5 Equation (9) gives the formula for the output layer: b) The error signal back propagation phase. The neuron error in the output layer is shown in (10): Equation (11) gives the formula for calculating the total network error E: Expanded: According to the BP principle, the amount of change in weight is proportional to the negative gradient of the error, which is obtained by the gradient descent algorithm: where The weight is updated to Above is the calculation process of the BP algorithm. Then, a set of optimal weights is trained to be applied to the control system.

c: BPNN CONTROL
When K 2 in Fig. 2 is switched on, the BPNN controller is applied, and the connection weights in the BPNN controller are trained through the above methods. Apply already trained weights to control other PD states (PD1 -PD3) which differ from the training data set. There are also two important VOLUME 8, 2020  parameters that affect network performance. They are the learning rate (η) and the number of neurons in the hidden layer (n hidden ), whose effects can be further explored. In this work, η and n hidden will be discussed in section III.

C. EVALUATION CRITERION 1) SYNCHRONIZATION INDEX (S I )
Experimental studies have shown that neurons in a normal state rarely fire synchronously, whereas neurons in a Parkinsonian state fire with a significant degree of synchronization [2]. Therefore, we define a measure of synchronization. When ten or more neurons simultaneously fire within a certain time window, this is defined as firing synchronously, and the average value over a period of time is defined as the synchronization index. Specifically, this work selects the firing activity during 4s to 5s to calculate the synchronization index. The quantity of synchronous occurrences is divided by the window length of 1s to give the synchronization index.

2) SUMMATION OF BETA BAND POWER (Q)
To evaluate the spectral characteristics of the different neural states, especially in the beta (12-35Hz) frequency band, this work also calculates the sum of power in the beta band. Based on the Fourier transform of the mean field signal of each nucleus (y I ), the power spectrum can be generated. The power for each frequency P for the nucleus I is calculated as follows: where f s is the sample frequency of the power spectrum, being set to 10 4 Hz here. t n represents the current time in the simulation process. Using the last 1s of the simulation results to calculate the power spectrum can ensure a reliable and stable result.
The total power in the beta band for nucleus I with its average value for the cortex-basal ganglia-thalamo are defined as follow: where Q I represents the sum of power in the beta band for nucleus I and N is the total number of neurons of the nucleus I .

3) ENERGY CONSUMED INDEX (E CI )
E CI is refer to the energy delivered to the tissue in this paper. The energy consumption is an important parameter to be considered. The energy consumption index is expressed by the following formula: where I DBS represents the stimulation current, and the total number of the neurons receiving the stimulation is N . Simulations are implemented in MATLAB 2018a (The Mathworks, Natick, MA, USA) with equations solved using the forward Euler method with a time step of dt = 0.1ms.

III. RESULTS
The initial membrane of each neuron was randomly distributed between [−80mV, −30mV], and the initial membrane recovery variable was randomly distributed between [-10mV, −1mV]. In order to remove the influence of the initial conditions on the simulation results, the simulation data of the first 1000ms was removed from the quantitative calculations, and all the results are averaged from 100 independent runs.
In order to reflect the robustness of the BPNN controller, three new configurations of PD states are selected which differ from the training data set: g Gpe→Gpe = 0.02 for PD1, g Gpe→Gpe = 1.2 for PD2, and g Gpe→Gpe = 0.6 for PD3.

A. MODEL VERIFICATION 1) FIRING RATES
The firing rate of each neuron was calculated for the time window from 1000ms to 5000ms, and then an average firing rate was calculated across all neurons in a given nucleus.
The changes in the firing rate of different nuclei in normal and three types of PD states are shown in Fig. 5. The firing rate of neurons in the STN and GPi nuclei increases, whereas that of the GPe decreased, when the state changes to the PD from normal. The simulation results of the model are consistent with experimental results [36]- [38]. Experimental studies of 6-OHDA-treated rats showed that the firing rate of STN and GPi nuclei increased during the transition from normal state to PD state [39], [40], but the firing rate of GPe neurons is reduced [41]. From the perspective of firing rate, the model reproduces the findings obtained in the animal model.

2) NEURAL FIRING SYNCHRONIZATION
A set of spatiotemporal firing pattern plots are presented in Fig. 6 to characterize the neural synchronization, where the ordinate represents a total of 100 neurons. We can see from Fig. 6, in the normal state, the synchrony of the neurons of each nucleus is low, whereas, in the Parkinsonian state, the phenomenon of synchronous firing occurs in all nuclei. Synchronization of the neurons in the PD state is significantly enhanced. This is consistent with the observation in the animal experiments [2].

3) NEURON POWER SPECTRAL DENSITY
By calculating the power spectral density of different nuclei in the normal and PD states, the established model in this work can be further verified [34]- [36]. Fig. 7 shows the average of the power spectral densities for each nucleus. By performing power spectral analysis, we can find that the peaks in the beta band  are enhanced in each PD state, as compared with the normal state.

B. PID CLOSED-LOOP DBS CONTROL RESULTS
In the PID controller, the controller parameters are set as K p = 10, K i = 0.01 and K d = 0 after many trials according to PD1 state. We can see from Fig. 8 that the PID controller can effectively alleviate the PD1 state. However, the synchronization phenomenon can also be observed in PD2 and PD3 states with control. Without changing the controller parameters, the PID controller can neither effectively control PD2 state nor PD3 states. Thus, it can be concluded that the PID controller cannot control different PD states under a fixed set of parameters. However, tuning of the PID parameters is undesirable because it is time consuming and full of challenges. Therefore, this paper hopes to find a method that can effectively control different PD states without changing the controller parameters.

C. BPNN CLOSED-LOOP DBS CONTROL RESULTS
The learning rate and the number of hidden layer neurons of the BPNN controller are significant parameters affecting controller performance, so it is important to choose their values reasonably. Their effects on synchronization and energy consumption are explored below, and reasonable values are determined. Fig. 9 explores the effect of n hidden and η on the S I . The n hidden is varied from 1 to 10, and the learning rate is varied from 0.01 to 0.05, respectively. We can see from Fig. 9 that when n hidden is 5 or 6, and η is between 0.02 and 0.04, S I takes the minimum value.
Next, Fig. 10 explores the effect of n hidden and η on the E CI . We can see from Fig. 10 that when n hidden is 5 or 6, and η is between 0.02 and 0.04, E CI takes the minimum value.
From the heat map of Fig. 9 and Fig. 10, the influence of n hidden and η on the control effect can be seen. It can obtain the good control performance when n hidden is 5 or 6, and η is between 0.02 and 0.04. Therefore, in order to obtain the best control effect, n hidden is chosen to be 5, and η is chosen to be 0.03. Moreover, We can see that the optimal BPNN parameters are suitable to the different PD states.
The raster and power spectral densities of PD1 -PD3 under the BPNN control are shown in Fig. 11. We can see from Fig. 11(a) to (c) that, the synchronizations of the neurons in each nucleus can be greatly reduced in PD1, PD2, and PD3 states under an identical BPNN controller. From the power spectral density plots, we can see that the BPNN controller effectively suppresses the beta rhythm of each nucleus in all PD states. Therefore, the obvious advantage of the BPNN control strategy is its robustness. Facing the multiple PD states, BPNN controller needs no more readjusted. Thus, it can be concluded that the BPNN controller has better control performance than the PID controller facing the control of Parkinsonian states.
In order to further evaluate the control effects, control signals from the BPNN controller are shown in Fig. 12.
It can be seen that, for different PD states, the control signals convergence within a short time.

D. CONTROL COMPARISON
In order to present the advantage of BPNN control strategy, more performance indexes are given to compare. It can be seen from Fig. 13(a) and (b) that both the PID controller and the BPNN controller have achieved effective control of the PD1 state. However, from Fig. 13(c) that the energy consumption of the BPNN controller is 58.26% lower than that of the PID controller. It can be seen from Fig. 13(a) and Fig. 13(b) that although the PID controller can effectively control PD1, it cannot effectively control PD2 and PD3, but the BPNN controller has achieved better control effects for PD2 and PD3, and the energy consumption is low. This shows that the BPNN controller proposed in this paper has better control effect than the PID controller on the control problem of Parkinsonian state.
Without changing the controller parameters, the BPNN controller can successfully modulate other different PD states and consume less energy, but the PID controller cannot. This shows that the BPNN controller has an advantage in the control problem of the Parkinsonian state.

IV. CONCLUSION
Based on the physiology of the brain and the firing mechanisms of neurons, a computational model of the basal ganglia-thalamo-cortical circuit is constructed. In view of the limitations of clinical open-loop DBS, this work constructs a closed-loop DBS system that can adjust the stimulus in real time according to the condition of a patient to achieve adaptive control. The research contents and results of this work are as follows: (1) Based on the physiological structure, the basal ganglia-thalamo-cortical computational model is constructed, and the normal and Parkinsonian states are simulated by adjusting model parameters. (2) Closed-loop DBS stimulation models are constructed to mitigate the Parkinsonian state. In the closed-loop DBS system, the BPNN controller can modulate different PD states without changing the controller parameters, but the PID-type controller cannot. Additionally, BPNN controller achieved better control performance by observing S I , Q and E CI . By comparing the energy expenditure of the two types of control strategies, it can be observed that the BPNN controller consumes less energy, which addresses a problem of excessive energy consumption in open-loop DBS systems and may prolong the life of the DBS battery, reducing the risk of surgery for frequent battery replacement.
However, The basal ganglia -thalamo -cortical computational model constructed in this work may not be very complete. In order to more realistically reflect the state of the neurons, it is expected that the striatum will also be constructed in the future. This work does not explore the basic mechanism of PD from an ionic channel level. In the future, the pathogenesis and treatment mechanism of PD will be further explored. How to mitigate the Parkinsonian state by the control strategy in finite-time maybe another interesting work. In the future, we will further explore to mitigate the Parkinsonian state in finite-time [42]- [44].
CHRIS FIETKIEWICZ received the B.S. degree in electrical engineering from the Messiah College, Mechanicsburg, PA, USA, and Temple University, Philadelphia, PA, USA, in 1991, through a joint program, and the Ph.D. degree in computer science from Case Western Reserve University, Cleveland, OH, USA, in 2010.
He was with the industry as an electrical engineer and a software engineer. He was a Postdoctoral Associate with The State University of New York. He uses techniques from computational neuroscience to better understand the nervous system. His current research interests include biological modeling and neuroinformatics. He is currently the Nord Professor of engineering with faculty appointments at the Department of Electrical Engineering and Computer Science, the Department of Mechanical Engineering, and the Department of Biomedical Engineering, Case Western Reserve University. His research interests include stability and control of nonlinear systems with applications to large-scale electric power systems, and nonlinear filtering with applications to monitoring, fault detection, diagnosis, and reconfigurable control. He is also a Fellow of the American Institute of Medical and Biological Engineers. VOLUME 8, 2020