Adaptive parameter modulation of deep brain stimulation in a computational model of basal ganglia–thalamic network

Deep brain stimulation (DBS) has proven to be an effective treatment for Parkinson’s disease (PD). Adaptive control strategies offer the potential to improve efficacy, limit side effects and save battery consumption via reducing the total amount of stimulation delivered. However, the mechanisms underlying the beneficial effects of DBS for PD remain poorly understood and are still under debate, which has hindered the development of closed-loop DBS. And during the chronically implanted phase, adaptive DBS needs to be further improved to maintain its advantages. In the design of new adaptive DBS, more insights into inaccuracies when establishing mathematical basal ganglia model, unknown external disturbance signal and dynamics of focal area should be considered. A controlled auto-regressive moving average is used as the representative description of stimulus–response relationship based on recursive extended least square method, where stimulation signal is applied to subthalamic nucleus (STN) and the feedback signal is selected as local field potential signal from internal segment of the globus pallidus (GPi). The generalized minimum variance algorithm is used for online update of stimulation frequency and amplitude in a closed-loop manner. Simulation results illustrated the efficiency of the proposed closed-loop stimulation methods in disrupting the aberrant beta-band (12–35 Hz) oscillation and restoring normal firing pattern when compared with the PD state. Robustness of the adaptive algorithm was further verified through dynamic change of illness state.


Introduction
Parkinson's disease (PD) is the second most common neurodegenerative disorder associated with the degeneration of dopaminergic neurons in the substantia nigra pars compacta [1][2][3][4][5][6]. Up till now, more than ten million people are suffering from the motor symptoms such as bradykinesia, rigidity, rest tremor and postural and gait impairment worldwide [7,8]. Conventional deep brain stimulation (DBS) has significantly improved the treatment of PD. Continuous highfrequency electrical stimulation is delivered to specific subcortical targets, including the subthalamic nucleus (STN) and the internal segment of the globus pallidus (GPi) [9][10][11][12][13]. However, the efficacy of the current open-loop DBS is restricted in terms of side effects, battery consumption and parameter programming, which presents demands for a highly automated and intelligent adaptive parameter modulation system [14][15][16][17][18][19].
In recent years, closed-loop stimulation utilizing patients' clinical states as a guide to alter stimulation parameters as necessary has attracted more attention [20][21][22][23][24][25][26][27]. It provides an effective approach to optimize clinical benefits while minimizing side effects and battery depletion. A critical step in the development of such systems is the identification of signal features or ''biomarkers'' which have the potential to quantify clinical states [28,29]. To date, local field potential (LFP) has become a commonly adopted feedback signal for closing the loop. On the one hand, LFP could be directly recorded from the contacts of the implanted DBS electrode, thereby avoiding the consumption of additional external sensors. On the other hand, it has been proved that LFP is closely related to synchronized bursting activity in Parkinsonian state [26,[30][31][32]. Little et al. adopted an ''on-off'' pattern to deliver the stimulation signals, in which the stimulator was switched on only when the detected LFP beta amplitude exceeded certain threshold [22]. Recordings reported that transient aggravation of disease was caused after the implantation of electrode for PD patients, but from a chronically implanted phase the above ''stun effect'' will disappear. The adaptive DBS outperformed conventional DBS in terms of efficacy and efficiency, becoming a good starting point for algorithmic control. However, it is worth noting that, prior to and during movement, synchronized beta oscillatory activity is reduced, which left the fixed LFP amplitude threshold not accurate for the representation of desired state. Thus, it is required to seek for solution to the problem by designing new adaptive DBS algorithm to modulate the stimulation parameters. Proportional control strategy has also shown potential in computational studies but was restricted due to its time-consuming parameter adjustment process [33]. However, it largely depended on the properties of controlled brain network, which left still a long way from clinical implementation. In practice, the precise mathematical model of the Parkinsonian focal area could not be established [34]. Besides, individual variation and illness progress need further to be characterized to enhance the effectiveness and efficacy of DBS [35][36][37]. The above issues put forward exacting requirements for the development of closed-loop strategy. To address the problem, a self-tuning controller, combining the recursive method of parameter estimation and closed-loop control algorithm, is proposed to build a real-time computer control system for the realization of adaptive adjustment of DBS parameters (including frequency and amplitude). Considering the existence of intrinsic time delay between excitatory or inhibitory neurotransmitter transmission within basal gangliathalamus (BG-Th) network, it is necessary to predict the control effect in advance to avoid suboptimal control effect [38][39][40]. Due to the characteristics of simple in algorithm and easy to implement, minimum variance control has become the foundation of selftuning control algorithms [41], which could be adopted in our work to develop a predictive closedloop system.
The rest of the paper was organized as follows. In Sect. 2, we introduced the models and methods. Firstly, we detailed the computational model of BG-Th network and how DBS signal was applied to neuron model. Next, we determined the selection of feedback signal and then generated input-output data for the identification of the stimulus-response relationship. Based on the above, we derived the real-time DBS stimulation signal and formed a consummate closed-loop control framework. Computer simulation results are presented in Sect. 3. Finally, the discussion and conclusion are given in Sect. 4.

Methods
In this section, we started by describing a previously proposed mathematical model of Parkinsonian focal zone, basal ganglia (BG). In this paper, it served as the data generator and stimulation testbed. Since Parkinsonian behaviors were related to rhythmic firing patterns with oscillations at beta band [29], the relationship between the applied stimulus signal and the calculated LFP output was established through real-time online identification. Based on the identified model, the closed-loop controller was developed and designed to select optimal stimulation to improve Parkinsonian state as shown in Fig. 1.

Model description
The biologically plausible mathematical BG-Th model presented by So et al. [42] was modified from classical Rubin-Terman model [43]. Compared with the original Rubin-Terman model, the change of ionic current types and the network topological structure prompted the firing properties better match experimental data. Therefore, it was an excellent platform to study effects of neural activation or silencing and to test the adaptive control algorithm in this paper.
The BG-Th network was formed by excitatory or inhibitory synaptic connections among four primary nuclei including STN, GPe (external segment of globus pallidus), GPi and thalamus (Th). Each STN neuron received inhibitory connections from two GPe neurons. Each GPe and GPi neuron received excitatory inputs from two STN neurons and inhibitory inputs from two GPe neurons. Each Th neuron received inhibitory input from one GPi neuron. All connections followed a structured and deterministic pattern and are illustrated in Fig. 2. The number of cells in each nucleus was 10, which was proven to be sufficient to generate same trends of normal and Parkinsonian states when compared with larger number of cells in each population. The single-compartment Hodgkin-Huxley-type equations were used to describe firing dynamics of STN, GPe, GPi and Th C m dv STN dt ¼ ÀI l À I Na À I K À I T À I Ca À I AHP where C m ¼ 1lF=cm 2 is the membrane capacitance, and v i ði 2 STN,GPe,GPiÞ represents the membrane potential. On the right side of above equations, Fig. 1 The diagram of CARMA model-based generalized minimum variance control system Fig. 2 The layout of BG-Th network and details of synaptic connections among STN, GPe, GPi and Th nucleus within the network, where red arrows denote excitatory connections, and black arrows denote inhibitory connections I l ; I Na ; I K ; I T ; I Ca ; I AHP were, respectively, the leak current, the sodium current, the potassium current, the low-threshold calcium current, the high-threshold calcium current and the after hyper-polarization potassium current. I app i ði 2 STN,GPe,GPiÞ represented constant positive bias current that inputs from unmodeled brain regions, such as striatum. The transition from normal state to Parkinsonian states was modeled by changing the value of I app i .pd ¼ 0 represented the normal state, and pd ¼ 1 represented the PD state, and thus, we got the following relationship Besides, I a!b ða; b 2 STN,GPe,GPiÞ describes the synaptic connections from pre-synaptic nucleus a to post-synaptic nucleus b, where g a!b and E a!b represent the maximal synaptic conductance and reversal potential, respectively, and P j s j a is the sum of synaptic variables. For synaptic projections that stem from GPe (from GPe to STN, from GPe to GPi, from GPe to GPe), the synaptic variables took the form as where H 1 ðvÞ ¼ 1=ðs þ expðÀðv þ 57Þ=2ÞÞ. In addition, for other synaptic projections (from STN to GPe, from STN to GPi, from GPi to Th), the dynamic of synaptic variables was calculated as follows where uðtÞ is an indicator of an action potential. That is, if presynaptic neuron produced an spike that crossed the threshold of -10 mV, uðtÞ ¼ 1. Otherwise, uðtÞ ¼ 0: In addition to ionic currents and synaptic currents, for Th neurons, I C represented a series of excitatory current pulses from cortical neurons and followed a gamma distribution with the average rate of 14 Hz and the variation of 0.2, the amplitude and pulse duration of which were set to 3.5 lA=cm 2 and 5 ms, respectively. Moreover, the implementation of DBS stimulus took the form as Conventional effective stimulus parameters were usually set to be A DBS ¼ 300lA=cm 2 , d DBS ¼ 0:3 ms and q DBS ¼ 7.69 ms (130 Hz), while in this paper, values of DBS parameters were modulated based on real-time calculation from closed-loop controller.
For the case of frequency modulation, For the case of amplitude modulation, where f max and f min are the maximum and minimum bounds of the modulated DBS frequency, A max and A min are the maximum and minimum bounds of the modulated DBS amplitude, and calculation of uðkÞ was introduced in the following part. Details of equations and parameters are given in Table 1, where membrane potentials, synaptic conductance, currents were expressed in mV, mS/cm 2 and lA=cm 2 , respectively. Simulations were implemented in MATLAB R2019a (The Math Works, Natick, MA), and the equations were solved using the forward Euler method with a time step of 0.01 ms.

Identification of the stimulus-response relationship
PD was a progressive neurodegenerative disease. Selection of appropriate biomarker to capture ongoing clinical changes was a critical step. It has been proved that beta-band oscillation of LFP signal was related to motor symptoms of bradykinesia and rigidity across patients. Recent research presented the evidence that LFP signal directly recorded from DBS electrodes appears to be a promising source of feedback signal. Based on the mathematical BG-Th network, the LFP signal of GPi nucleus [44] could be computed by which was an approximation of the sum of transmembrane currents. However, due to the characteristic of highly nonlinearity within BG-Th network, it presented huge challenges in the design of closed-loop control algorithm. A recursive extended least square (RELS) method was proposed to realize online real-time parameter estimation, and therein, the model structure could be delineated by the controlled auto-regressive moving average (CARMA) model where uðkÞ represents the input DBS stimulus input (to STN), yðkÞ represents the output LFP signal (from GPi) at time kT s (T s = 0.01 ms is the sampling interval), nðkÞ is Gaussian white noise series representing the sum of the stochastic perturbation and unmodeled dynamics, and d is the time delay. Coefficients of polynomial ð11Þ remained to be determined by detectable input and recorded output data. Firstly, the CARMA model was written in least squares form where uðkÞ ¼ ½Àyðk À 1Þ; :::; Àyðk À n a Þ; uðk À dÞ; :::; uðk À d À n b Þ; nðk À 1Þ; :::; nðk À n c Þ T h ¼ ½a 1 ; :::; a n a ; b 0 ; :::; b n b ; c 1 ; :::; c n c : Due to the uncertainty of nðkÞ in uðkÞ, it was substituted using the estimation value n _ ðkÞ ¼ yðkÞ À y _ ðkÞ ¼ yðkÞ À u _ T ðkÞh. Then, least squares estimation was performed to minimize the performance evaluation The root-mean-square error was calculated to quantify identification precision between actual value and estimated value. In consideration of model complexity and computational efficiency, Akaike information criterion (AIC) served as a reference to help us find a relative optimal fitting model within limited order range, where K represents the total number of parameters to be estimated, N is the sample capacity, and L ¼ ÀN=2 Ã ðlnð2pÞ þ 1 þ lnðe 2 =NÞÞ.

Controller design
To automatically adjust DBS parameters, a real-time computer control system was established by combing the recursive method of parameter estimation and selftuning control method. Considering the existence of time lag between DBS action and LFP output, the generalized minimum variance control (GMVC) predicted the output in advance and then designed the control law to guarantee the stability of BG-Th network [45].

> : ð24Þ
By this time, the cost function was expressed in The first term was omitted due to its independence with past input and output, and then, we took the partial derivative of rest terms in Eq. (25) with respect to uðkÞ 2 Py Ã ðk þ dÞ À Ry r ðk þ dÞ ½ oPy Ã ðk þ dÞ ouðkÞ Substituting oPy Ã ðkþdÞ ouðkÞ ¼ b 0 , oQuðkÞ ouðkÞ ¼ q 0 into Eq. (26), the generalized minimum variance control law was obtained as uðkÞ ¼ Ry r ðk þ dÞ À GPyðkÞ where the control signal at time kT s is obtained and applied to the BG-Th model target.

Simulated BG-Th network firing pattern
The transition from normal state to Parkinsonian state was achieved by changing the level of bias currents to STN, GPe and GPi neurons. As shown in Fig. 3a, b and c, the STN, GPe and GPi neurons exhibited regular spiking in the normal physiological condition, while in the PD state, a rhythmic bursting fashion occurred, especially in the GPe and GPi neurons as shown in Fig. 3d, e and f. Then, we measured the relay ability of Th neurons in the reference of cortex input. Error index was a valid proxy to quantify the reliability of Th neurons. A spike of Th neuron was defined when the transmembrane voltage crossed the threshold of -40 mV. Correspondingly, an error was defined when a spike failure occurred in response to cortex input. N error represented the total number of errors (including missing spikes, rebound bursts and spurious events) in thalamic transmission, and N c summed the number of cortex input pulses. As shown in Fig. 4, more firing errors occurred due to rhythmic pattern of GPi neurons in the pathological condition. Statistical results revealed a significant increase of the error index EI from normal condition (EI normal ¼0:0093 AE 0.0011) to PD condition (EI PD ¼0:4865 AE 0.0787), as shown by the first two bars in Fig. 5. Since it has been proved in clinical that  [47,48] but stimulation with frequencies below 50 Hz might even exacerbate symptoms [49][50][51][52], we explored the frequency-dependent effectiveness of DBS through the way of quantifying the reliability of Th. With the increase in frequency of DBS, the trend of error index in yellow bars provided sound evidence of the above facts. Besides, the LFP signals of GPi under normal condition and PD condition are presented in Fig. 6a and b, respectively. As shown in Fig. 6c, power spectral analysis manifested an obviously beta-band peak under PD condition when compared with normal condition.

Performance of stimulus-response identification
A CARMA model was identified to establish the relationship between DBS stimulation signal and the LFP response signal. Based on the obtained input and output data from the BG-Th network, it was necessary to select appropriate order parameters by quantifying the estimation accuracy. The root-mean-square error between actual value and estimated value of LFP signal in GPi was calculated and compared among different orders of n a and n b . The trend of root-meansquare error and the AIC with the change of n a and n b is plotted in Fig. 7a and b, respectively. It was obvious that the fitting precision was improved with the increase of n a and n b , while AIC showed an obviously opposite trend. To reach a compromise between precision and efficiency, n a and n b were both set to 4. Parameters c 1 ; c 2 ; :::; c n c ½ were used to measure the influence of uncertain dynamics, and n c was set to 3 due to its minimal impact on both the indexes above.
In this case, evolution of estimation results is illustrated in Fig. 8a-c, where all identified parameters rapidly converged to stable values. Correspondingly, the identification performance is exhibited in Fig. 8d, where the red solid line displayed the original data recorded from BG-Th network output and the blue dotted line represented the estimated data generated from CARMA model. The result suggested that the CARMA model provided a good fit to the obtained input-output data.

Efficacy of GMVC strategy
Based on the identified parameters of the CARMA model, the frequency or amplitude parameter was calculated via the adaptive generalized minimum variance controller to complete the design of closedloop DBS system. As stated in Sect. 2, we manipulated the STN neurons as stimulation targets and obtained real-time feedback from the GPi neurons. For the case when closed-loop frequency stimulation was implemented (since t = 1000 ms) to steer DBS waveform, rhythmic bursting firing pattern of GPi was obviously suppressed, which was illustrated by comparison of LFP waveform signal before and after stimulation in Fig. 9a. Figure 9b indicates the change of stimulation frequency by the sparse and intensive distribution stimulus pulse, where f max and f min are set to 200 Hz and 50 Hz to simulate clinically effective frequency bands. As shown in Fig. 9c, quantitative analysis through power spectrum indicated the disrupted pathological frequency structure after closed-loop DBS amplitude modulation. Similarly, for the case when DBS was applied to STN neurons with real-time amplitude modulation (f ¼ 130 Hz), suppression of rhythmic oscillated firing pattern and disruption of pathological spectrum structure were realized since t = 1000 ms as shown in Fig. 10. Further analysis of closed-loop performance was developed in terms of interspike interval (ISI) changes. As shown in Fig. 11, abscissa presented ISI distribution and ordinate empirical cumulative distribution function. Regular spike firing in normal condition was always accompanied by uniform distribution of small ISI, while burst firing pattern was accompanied by dispersion of ISI. Specifically, as shown in Fig. 11, ISI distributed mainly around 20 ms in normal state, but in PD state the distribution of ISI was focused around 10 ms and 60 ms, representing the strong polarization induced by increase in pathological burst firing [46]. Under both closed-loop frequency modulation and amplitude modulation, the pathological condition of ISI distribution was gradually recovered with the tendency toward normal distribution.
As described ahead, the Parkinsonian state was simulated by decreasing the constant bias currents applied to the STN, GPe and GPi neurons. Directly switching I app i from a higher value to a lower value apparently ignored the possible intermediate states resulting from diverse courses of disease caused by individual differences, and thus, we defined a dynamic parameter pd to quantify disease severity. Correspondingly, the robustness and adaptability of GMVC strategy on closed-loop control performance were tested by continuously changing pd from 0 to 1, wherein the Beta index was defined as the ratio of integrated spectral power within 12-35 Hz to integrated spectral power within 0-50 Hz, to diagnose the existence of beta peak and to distinguish normal state from PD state. As shown in Fig. 12, the red error bar presented the trend of Beta index with the increase of pd. Strategy of real-time closed-loop frequency modulation and amplitude modulation was implemented when pd increased from 0.1 to 1. In both cases, the increased Beta index was almost pulled back to normal state, denoted by gray error bar with yellow round markers and gray error bar with green square markers, indicating the restoration of the spectral features.

Discussion and conclusion
The objective of this study was to design a closed-loop control strategy to suppress Parkinsonian state. To advance the progress of neurophysiological experiments and to reveal the principles of information processing together with various functions of the nervous system, we focused on the model research to inspire new treatments for neurological disorders. In this paper, the computational BG-Th network consisting of STN, GPe, GPi and Th neurons was constructed according to Hodgkin-Huxley-type equation, with parameter settings and synaptic connections in conformity to anatomy and physiological properties.
Since the findings are that the rhythmically oscillatory activity of LFP is a remarkable PD-related feature, LFP of GPi neurons was adopted as a feedback signal that distinguishes normal condition from PD condition. Due to the lack of adequate physiological details plus the uncertainties in noise or disturbances, it was utilized to adopt a black-box CARMA model to remedy existing limitations. Therefore, relationship between stimulation signal and neuronal response was expected to be described through the identification method. Further, adaptive reprogramming DBS stimulation parameters in reducing the adverse effects of continuous high-frequency DBS were tested through the proposed self-tuning GMVC control strategy. Numerical simulation results illustrated the possibility of mitigating intractable Parkinsonism by GMVC, including the suppression of abnormally enhanced beta-band oscillation and the regulation of interspike interval distribution. Firstly, the intrinsic transmission time delay between STN and GPi occupied an important and non-negligible position, which was considered in the design of the closed-loop system. Secondly, it was hypothesized to represent possible dynamics of disease progress by adjusting the parameter pd from 0 to 1, and simulation results indicated that the GMVC strategy reached remarkably excellent controlling performance. Although robustness of the proposed adaptive control algorithm was proved, the modified computational model adapted from Rubin-Terman model was still far from enough to reveal the real brain network and it inevitably ignored considerable physiological knowledge. Therefore, the potential application of GMVC in future clinical DBS requires intensive consideration on practical issues, including which is the most appropriate biomarker for the feedback variable and how biomarkers vary across patients and progress over c Comparison of the power spectral analysis before and after the external stimulation. Controller parameters were set the same as the case of frequency modulation