A biophysically-based model of the cortex-basal ganglia-thalamus network (Kumaravelu et al., 2016), modified from the original Rubin-Terman model (Rubin and Terman, 2004), was adopted as a platform to develop and evaluate the controllers. The model included representations of neurons in cortex (CTX), striatum (STR) [sum of direct striatum (dSTR) and indirect striatum (idSTR)], STN, globus pallidus (GP) (sum of externa part GPe and interna part GPi) and thalamus (TH), and all the nuclei were interconnected through either excitatory or inhibitory synaptic connections to form a network. Each nucleus contained 10 single-compartment model neurons (Figure 1).

The neurons of the STN, GP, TH, and STR were modeled using Hodgkin-Huxley (HH) type equations:

where C represented the membrane capacitance and was set to 1μF/cm² for all cell models, v_i(i ∈ {STN,GP,TH,STR}) represented the transmembrane potential of the corresponding model neuron and was expressed in mV. I_Na,I_K,andI_l were the sodium current, potassium current and non-specific leak current, I_L,I_T,andI_Ca were, respectively, the L-type, T-type, and high-threshold calcium current, and I_a,I_m,I_Cak,andI_AHP were, respectively, A-type, outward M-type potassium, calcium-dependent and after-threshold potassium current. Parameters and equations of the ionic currents are provided in Table 1. Therein, a, b, c, d1, d2, h, m, n, p, q, and r were activation or inactivation variables, and the gating kinetics took the form

where X represented one of a, b, c, d1, d2, h, m, n, p, q, or r. Steady-state gating variables were calculated using

where w_X and σ_X were the half voltage and slope, respectively. Gating kinetics for STR took the form

I_syn represented the sum of synaptic currents, with each projection from presynaptic neuron α to postsynaptic neuron β(α,β ∈ {CTX,STR,STN,GP,TH}) given by I_α→β = g_α,β×(v_β−E_syn)×S, where g_α,β described the maximal synaptic conductance, and E_syn represented the reversal potential (uniformly set as −85 mV). An alpha synapse S was used to model the synaptic dynamics

where t_d was the synaptic transmission delay, and τrepresented the time constant of 5 ms. Further, bias currents I_{app_i}(i ∈ {GPe,GPi,TH}) represented other synaptic inputs that were not described explicitly in this model.

The dynamics of the CTX neurons were described based on the model developed by Izhikevich (2003)

where v_CTX represented the transmembrane potential and u_CTX was the recovery variable. The time scale of the recovery variable u_CTX was chosen as a = 0.02 and a = 0.1, respectively, for excitatory-CTX (eCTX) and inhibitory-CTX (iCTX) neurons. If the transmembrane potential of CTX neuron exceeded 30 mV, then v_CTX was set to a resting potential equal to -65 mV and u_CTX was set to u_CTX + d (for eCTX, d = 8; for iCTX, d = 2). I_interCTX represented the reciprocal synaptic current from iCTX neurons to eCTX neurons or that from eCTX neurons to iCTX neurons, and I_TH→CTX was the synaptic input from TH.

The Parkinsonian state was simulated by the adjustment of model parameters implemented using a parkinsonism variable, pd, where pd = 0 and pd = 1 were defined as the healthy state and full Parkinsonian state, respectively. The M-type potassium current in striatal neurons was reducedg_m = 2.6−0.9×pd, cortico-striatal coupling strength was decreased g_CTX,STR = 0.07−0.044×pd and coupling strength between GPe neurons was increasedg_GPe,GPe = 0.0125 + 0.0375×pd. To quantify the difference between the healthy and Parkinsonian states, changes of firing rates and firing patterns of model neurons were analyzed. It was considered that one neuron produced a spike or action potential once its transmembrane potential was greater than the threshold V_thre = −20mV, with the time of crossing the threshold defined as the firing time. The average firing rates were calculated based upon the firing time during the entire simulation period. In addition, spike synchrony that characterizes the dynamic patterns in each population of model neurons was measured. Defining vⁱ(t) as the membrane potential time course of the jth neuron from a population of n neurons, then we could average over the population V⁢(t)=1n⁢∑i=1nvi⁢(t). The variance of membrane potential vⁱ(t) and the variance of the averaged membrane potential V(t) were expressed as σvi2=⟨[vi⁢(t)]2⟩t-[⟨vi⁢(t)⟩t]2 and σV2=⟨[V⁢(t)]2⟩t-[⟨V⁢(t)⟩t]2, respectively (⟨⋯⟩t=1T∫0T⋯dt referred to the average value of the variables within the time of T), and the level of synchrony was calculated according to the following equation,

where χ was normalized between 0 and 1, with χ = 0 indicating neurons within one population fire out of sync and χ = 1 indicating neurons within that population discharge synchronously.

Parkinson’s disease is characterized by diverse changes in neuronal activity, and single neuron action potentials, electrocorticograms, LFPs, and electroencephalograms have been considered as feedback control signals for closed-loop DBS (Hoang et al., 2017). The LFP generated by model GPi neurons was adopted as the feedback signal for closed loop control. A simple average of transmembrane potentials was adopted to calculate the LFP of the modeled population due to its ability to capture subthreshold activity and thereby reflect oscillatory phenomena (Pettersen et al., 2012; Mazzoni et al., 2015). Expression of the GPi LFP was given as

Here, n = 10 represented the total number of GPi neurons and v_{GPi_i} corresponded to the transmembrane potential of the ith GPi neuron. The power within particular frequency bands of the LFP signal was determined from power spectra using the Chronux neural signal analysis package [length of moving window 1 s, step size 0.1 s and tapers in the form of [3 5] (3 is the time-bandwidth product and 5 is the number of tapers to be used)], and the beta band power was defined as the total power over 13–30 Hz.

Our goal was to design an adaptive closed-loop controller to adjust automatically the stimulation signals delivered to the model neurons of the GPi based on the beta LFP activity calculated from the model GPi neurons as the feedback control signal. We defined the stimulation signal as I_sti and delivered it directly to each GPi neuron, and the resulting transmembrane potential was expressed as

I_sti was constructed by using the controller output u(t) to construct a variable frequency pulse train stimulation signal [amplitude of 300μA/cm², pulse duration of 0.3 ms and period of 1,000/u(t) ms (frequency of u(t))]. After the end of each stimulation period, we recalculated the beta power of the LFP signal and repeated the above steps to update continuously the optimal controller output. The transmembrane potentials of model neurons in that other nuclei were unaffected by direct electrical stimulation but their activity was influenced during stimulation via either excitatory or inhibitory synaptic connections.

The classical error-based PID control law has the form

where e(t) = y_d(t)−y(t)(y_d(t) and y(t) represented the beta power in LFP signals of the GPi from healthy control and controlled Parkinsonian states, respectively). The performance of the PID controller depends greatly on selecting the appropriate gains, and this can be a time-consuming manual process. As the LFP was a highly dynamic variable subject to large changes, differential action might amplify noise interference. Therefore, the differential term was omitted, and both proportional (P) and proportional-integral (PI) controllers were designed to minimize the error between the desired and measured beta band power in the LFP.

A stable self-tuning controller was designed using a dynamic RBF network. Figure 2 is a block diagram of the RBF supervisory control system and a schematic diagram of the RBF network. The neural network-based controller acted as a feedforward controller, by building an inverse model of the controlled plant. The input layer, hidden layer and output layer determined the structure of the RBF neural network. x = [x₁,x₂,…,x_i,…,x_n]^T (n was the number of input layer nodes) represented the network input vector and h = [h₁,h₂,…,h_j,…,h_m]^T (m was the number of hidden layer nodes) represented the hidden layer output. Each hidden layer node had a central value c_j, the Euclidean distance of which to network input x_i was described as ||x_i−c_j||. As well, each hidden layer node was an arithmetic element with activation function given by

The center vector c = [c₁,c₂,…,c_m] and the width vector b = [b₁,b₂,…,b_m]^T determined the influence of the Gaussian function, where the width of the Gaussian basis function directly influenced the mapping capability to network input, while the center value correlated with its sensitivity to network input.

If we set the weight vector as w = [w₁,w₂,…,w_m]^T, then the network output can be obtained as

The structure of the RBF network was selected as 1−n_m−1, that is, containing 1 input layer node, n_m hidden layer nodes and 1 output layer node. The number of hidden layer nodes was set to 11, network weights were initially set to random values between 0 and 1, and parameters of Gaussian function were set as c = [−2,−1,0,1,2]^T, b = [5,5,5,5,5]^T. The weights of the RBF controller were continuously adjusted on-line to make the feedback error e(t) approach zero, which equated to u_p(t) approaching zero. Consequently, the RBF controller gradually occupied the leading position and even replaced the function of the P/PI controller. The RBF network error index was designed in the form of E⁢(t)=12⁢(up⁢(t))2, to lead u_p(t)(e(t)) to converge to 0. Considering that the total controller output was the sum of the traditional P/PI controller and the adaptive RBF controller u(t) = u_rbf(t) + u_p(t), the error index can be written as

According to the gradient descent method, the network weights were adjusted as follows

Further, applying the gradient descent method to the adjustment of c and b will optimize effective learning by the RBF network, thus we had

where η ∈ (0,1) represented the learning rate η = 0.30and α ∈ (0,1) represented the momentum factor α = 0.05.

A quantitative index of the control effect was defined as the root mean square error between the controlled output and the reference signal,

where N represented the sampling point of the feedback signal.

The transmembrane potentials of model GPe, GPi, and STN neurons in the cortical-Basal-thalamus network model are displayed in Figure 3. The Parkinsonian condition resulted in changes in both the rate and pattern of model neuron activity. As a result of excitation via the indirect pathway and hyper-direct pathway together with inhibition via the direct pathway, firing rates in the Parkinsonian condition increased in STN and GPi model neurons and decreased in GPe model neurons (Figure 4A), consistent with previous experimental studies (Kita and Kita, 2011). Moreover, increased spike synchrony was observed across all nuclei in the Parkinsonian state compared to the healthy state (Figure 4B). These indexes were analyzed using one-way analysis of variance (ANOVA), which revealed a significant difference in firing rate (STN: F = 901 and p < 0.001, GPe: F = 48 and p < 0.001, and GPi: F = 184 and p < 0.001) and synchrony index (STN: F = 239 and p < 0.001, GPe: F = 1554 and p < 0.001, and GPi: F = 62 and p < 0.001) between healthy state and Parkinsonian state. The model thus exhibited features of the pathophysiological neural activity occurring in PD and was a suitable testbed to develop and analyze closed-loop control strategies.

The GPi, which is clinically accessible for both recording of LFPs and delivery of DBS, was selected as the source of the feedback signal for closed loop DBS (Figures 5A,B). Compared to the healthy state, where the LFP exhibited little power in the beta band, the LFP in the Parkinsonian condition exhibited oscillatory activity, generating a significant peak in the power spectrum (Figure 5C). The LFP signals were filtered within the beta band to extract differences between healthy and Parkinsonian states, and the filtered GPi LFP activity in the healthy state (Figure 5D) served as a reference to guide the modulation of Parkinsonian neural activity by DBS, thus constituting a closed-loop system to suppress exaggerated beta oscillatory activity (Figure 5E).

We first quantified the relationship between the stimulation frequency and the beta power (Figure 6), where f = 0 is equivalent to the Parkinsonian state without DBS. Low frequency stimulation actually increased beta power, and the beta power was progressively reduced as the stimulation frequency was increased higher than 50 Hz. A beta power of 120 dB was calculated from the healthy state and set as the desired state, and the range of stimulation frequencies was between 5 and 200 Hz. The modulation of DBS was reformulated as a train of monophasic pulses with an initial stimulation frequency of 5 Hz, and the updated stimulation frequency was generated based upon the error between the actual and target beta power in the LFP from the GPi model neurons.

The performance of the P controller and PI controller in regulating the beta band LFP in the GPi with GPi DBS were compared (Figures 7, 8), where the blue dotted lines represented the healthy state and the orange traces represented the controlled PD state. Small k_p appeared to make no difference to the suppression of beta band activity, while large k_p caused strong oscillations between effective and ineffective suppression of beta band activity (Figure 7). The performance following addition of the integral controller with different combinations of k_p and k_i was assessed (Figure 8). The control of beta power in the GPi LFP was strongly dependent on the selection of k_p and k_i, where effective values (k_p = 0.5,k_i = 0.5) promoted the suppression of high beta power while ineffective values (k_p = 0.1,k_i = 0.01) did not contribute to improvement of the PD state. The search of proportional and integral gains through trial-and-error produced fluctuations in performance and created uncertainty about the effectiveness of closed-loop DBS, especially in the face of changes in the properties of the plant.

A substantial improvement in control performance was achieved by the supervisory algorithm, where suppression of the exaggerated beta power present in the Parkinsonian condition was achieved within 1 s (Figure 9A). The weights of the RBF network were adjusted in real-time, in response to the update of the beta power of the GPi (Figure 9B). As the beta power was gradually suppressed, the RBF network took over the leading role that the P/PI controller played in the initial control stage (Figure 9C). Ultimately, the DBS pulse repetition frequency was calculated as shown in Figure 9D, and the DBS signal is depicted in Figure 9E. Changes in beta power between the healthy state, open-loop 130 Hz DBS, P controller, PI controller and adaptive DBS were compared by calculating the root mean square error between the controlled output signal and the reference signal during the last 4 s of simulation. The averaged RMSE were 18.64, 112.41, 39.38, and 23.89, respectively, In cases of open-loop 130 Hz DBS, P control (k_p = 0.1), PI control (k_p = 0.5,k_i = 0.5) and combined P and RBF control, indicating that the improved supervisory algorithm drive the beta power to the target setting with a higher accuracy.

The closed-loop RBF network-based supervisory control algorithm achieved effective tracking of the target beta power, but several challenges required further consideration. First, Parkinson’s disease is a chronic and progressive disease in which the patient’s condition gradually worsens over time, and individual variation should be considered in model-based evaluations. Second, the reference waveform was the LFP signal from the healthy network model. Although this signal carries abundant information, it does not represent the variety of different disease states, for example during cycling of medication, and this potentially limits the generalizability of the control system.

To evaluate further the performance of the proposed closed-loop algorithm in the face of changes in the Parkinsonian state (Figure 10), the parkinsonism variable pd was randomly generated from 0 to 1. Controlled beta power gradually converged to the desired healthy signal after 1 s, demonstrating the adaptive capability of the RBF network across disease states. In addition, beta power exhibits dynamic changes, especially prior to and during movement, and thus tracking of time-varying beta power may be required to promote desired movement behavior. In the face of a time varying reference beta power signal switching at 1 Hz (Figure 11), the controlled beta power still followed the dynamic reference signal, albeit with substantial overshoot.