The computational model consists of four components: (i) neural activities, spikes, generated by the Vim network model in response to different DBS frequencies; (ii) motor cortex neural activities influenced by propagation of Vim–DBS effects to the motor cortex; (iii) spinal motoneuron activities impacted by neurons in the motor cortex; (iv) motor unit action potentials in the muscle fibers innervated by the spinal motoneurons.

Our computational model simulates the EMG signal in response to different DBS frequencies. Simulated EMG signals are used to calculate the feedback biomarker which controls the DBS frequency. The feedback control consists of three main parts: (1) biomarker identification, (2) computation of a system output from the biomarker, and (3) a closed-loop controller that implements the system output to update the DBS frequency.

The EMG simulation of our computational model is slow to implement: it takes more than 30 min when using MATLAB R2022b with a personal computer to simulate 10 s of the EMG signal. Thus, to facilitate the implementation speed of the model in a closed-loop control system, we need rapid EMG estimation to replace the direct EMG simulation. We implemented a polynomial method to estimate the EMG from the Vim firing rate, and the direct model-simulated EMG is used as the reference (i.e., reference EMG) for estimation. We implemented the MATLAB custom function polyfit for polynomial estimation, which gives a least-square fit of the polynomial coefficients. The estimated EMG is formulated as follows: y ^ t = ψ ζ x t + φ 0 + ε t ψ ζ x t = ∑ n = 1 N φ n ζ x t n  and  ζ x t = x t   −   x ¯ t s d x t , (6) where x t is the Vim firing rate simulated from our previous Vim network model (Tian et al., 2023b) and y ^ t is the estimated EMG. ζ x t is the standardization of x t ; x ¯ t and s d x t are the mean and standard deviation of x t , respectively. Then, ζ x t has a mean of 0 and a standard deviation of 1. The polynomial order is 25, and φ 0 , φ 1 , … , φ N , N = 25 are the polynomial coefficients. The R 2 statistic (Colin Cameron and Windmeijer, 1997) of the fit generally increases with increase in the polynomial order (Supplementary Figure S3), and 25 is a minimal order of the polynomial when R 2 converges. ε t is the Gaussian white noise with a standard deviation of 0.038 mV. We fitted the consistent polynomial coefficients (Supplementary Table S2) across data from different DBS frequencies, including 10, 50, 80, 100, 120, 130, 140, 160, and 200 Hz. Our previous work showed that consistent model parameters fitted based on concatenated DBS frequencies in a certain range—in this case, [10–200] Hz—can be consistently applied to unobserved frequencies (e.g., 25 and 180 Hz) in the same range (Tian et al., 2023a). Thus, we implement the polynomial-estimated EMG as the biomarker to control the DBS frequencies in the range [10–200] Hz.

The power spectral density (PSD) of the estimated EMG y ^ t (Eq. 6) is the system output for updating the frequency of the input DBS (10–200 Hz). In the system output, we consider PSD in the frequency band [2–200] Hz, which includes the frequencies of both DBS and EMG activities (Halliday et al., 2000; Hess and Pullman, 2012; Herron et al., 2017; Milosevic et al., 2021). The power density is calculated by the following equation (Miller, 2019; Liu et al., 2021), with a sampling frequency of 0.5 kHz: p f , t , u = 1 w ∫ − w 2 w 2 y ^ u t + s H s e i 2 π f s d s 2 ;  furthermore , (7)

H s = 1 2 1 + ⁡ cos 2 π s w .

where f represents the frequency of the EMG activities. H s is the Hann windowing function (Miller, 2019; Porat, 1997), and w is the width of the window, chosen to be w = 1 s to ensure stability (Liu et al., 2021; Miller, 2019). y ^ u t represents the estimated EMG in response to DBS with stimulation frequency u . The total power of EMG activities over the initial period T of the estimated EMG signal is then given as P f , T , u = ∫ 0 T p f , t , u d t . (8)

The EMG power thus computed can be used to analyze the EMG activities in the frequency domain in response to different frequencies of DBS.

In response to DBS frequency u , the PSD of EMG within the frequency band of [2–200] Hz is approximated as the sum of the power density ( p f , t , u , Eq. 7) from each integer DBS frequency (2, 3, 4, …, 200 Hz) in [2–200] Hz, and the system output z u is calculated as the approximated PSD in the initial T = 5 s of the estimated EMG: z u = 1 T ∫ 0 T ∑ f = 2 200 p f , t , u d t . (9)

Finally, a closed-loop controller updates the DBS frequency u so that z u is close to a specified control target β z (see the next section).

The simulation study for the Vim–DBS control system is schematized in Figure 2. The effects of Vim–DBS included direct axon activation and DBS-induced Vim firings, which were simulated from our previous Vim network model established based on clinical Vim–DBS data (Tian et al., 2023b). The Vim–DBS effects are propagated to the M1 neuron, which projects to the motoneuron in the spinal cord. The firings of the spinal motoneurons innervate the corresponding motor units in the muscle fibers and induce motor unit action potentials (MUAPs). The simulated EMG consists of a linear summation of MUAPs and a low-level Gaussian white noise. In the feedback control, in order to facilitate the implementation speed, we estimated the simulated EMG by a polynomial fit. Then, we computed the mean power spectral density (PSD) of such polynomial-estimated EMG as the system output. Finally, a proportional–integral–derivative (PID) controller updated the DBS frequency that brought the system output close to a specified control target.

DBS pulses are delivered to the neurons in the thalamic ventral intermediate nucleus (Vim) and also activate the axons projecting to neurons in the primary motor cortex (M1). The firing rate of the Vim neurons impacted by DBS was obtained from our previous Vim network model that reproduced experimental DBS data (Tian et al., 2023b). Spike trains are generated from the modeled Vim firing rate ( x t ) as a Poisson process and are propagated to the M1 neurons. The spikes from the M1 neurons are then propagated to the motoneurons in the spinal cord. The spikes from these motoneurons innervate the corresponding motor units in the muscle fibers and induce motor unit action potentials (MUAPs). The simulated EMG ( y t ) is a linear summation of MUAPs with additive low-level Gaussian white noise. In the feedback control, we estimate the simulated EMG by fitting a polynomial function ψ (Eq. 6), and this estimated EMG ( y ^ t ) is the biomarker used in the control. φ 0 is the constant term of ψ , ε t is the Gaussian white noise, and ζ x t is the standardization of the modeled Vim firing rate x t (Eq. 6). The system output z u is calculated as the mean power spectral density (PSD) over the initial T of the biomarker y ^ t (Eq. 9). Finally, a proportional–integral–derivative (PID) controller (Eq. 10) updates the DBS frequency u that reduces the error ( e ) between the system output z u and a specified control target β z .

Our model of the propagation of the Vim–DBS effects to the cortical M1 neurons was validated by using a 10-s block of neural spike data from single-unit recordings in M1 during 130-Hz VPLo-DBS in non-human primates, reported in Agnesi et al. (2015). The thalamic VPLo nucleus (ventralis posterior lateralis pars oralis, in the Olszewski atlas) in non-human primates is homologous to the Vim in humans (Molnar et al., 2005; Xiao et al., 2018). To assess the response of Vim to 130-Hz DBS, a peristimulus time histogram (PSTH) was calculated based on spike times occurring between 0 and 7.7 ms after each DBS pulse, with the PSTH further smoothed by an optimal Gaussian kernel (Shimazaki and Shinomoto, 2007; Shimazaki and Shinomoto, 2010) (Figure 3), both for empirical and simulated data. The standard deviation of the Gaussian kernel was 0.2 ms, which was obtained from optimization of the Gaussian kernel to best characterize the spikes using a Poisson process (Shimazaki and Shinomoto, 2007; Shimazaki and Shinomoto, 2010; Tian et al., 2023b). Our Vim–M1 propagation model’s generated M1 spike activity behaved similarly to that of the empirically recorded non-human primate M1, as indicated by the spike raster plot and PSTH firing rate analysis (Figure 3). We computed the R 2 statistic (Colin Cameron and Windmeijer, 1997) to compare model-simulated and experimental PSTHs, and R 2 = 0.728 represents a good model fit (Figure 3B). Thus, our Vim–M1 propagation model could reflect the M1 dynamics during 130-Hz Vim–DBS.

We simulated the EMG from our model, both during DBS-OFF and in response to Vim-DBS of different stimulation frequencies in [10–200] Hz (Figure 4).

The EMG simulation with DBS-OFF presented a typical tremor band (∼6 Hz) in the clinical EMG signals recorded from ET patients (Halliday et al., 2000; Hess and Pullman, 2012; Herron et al., 2017). During low-frequency ( ≤ 50 Hz) DBS, the amplitude of the simulated EMG is similar to (or slightly higher than) the DBS-OFF situation (Figure 4). Such a simulation is consistent with the clinical observations that low-frequency Vim-DBS ( ≤ 50 Hz) is often ineffective and can exacerbate the tremor (Ushe et al., 2004; Earhart et al., 2007; Pedrosa et al., 2013). The amplitude of the simulated EMG is lower than that of the DBS-OFF situation when the DBS frequency is ≥ 80 Hz (Figure 4). During high-frequency ( ≥ 100 Hz) DBS, the amplitude of the simulated EMG is clearly depressed compared with that in the DBS-OFF situation (Figure 4). Such a simulation is consistent with the clinical observations that high-frequency ( ≥ 100 Hz) Vim–DBS can worsen the tremor (Earhart et al., 2007; Ushe et al., 2004; Vaillancourt et al., 2003). The simulated EMG is mostly suppressed when the DBS frequency is ≥ 130 Hz (Figure 4). This is consistent with the fact that the stimulation frequency of clinical Vim–DBS is usually chosen to be ≥ 130 Hz (Ondo et al., 1998; Dowsey-Limousin, 2002; Dembek et al., 2020). We observed a short transient with a large amplitude in the simulated EMG during Vim–DBS, and the tremor intensity might be higher in the initial ∼200 ms after Vim–DBS onset (Milosevic et al., 2018; Yamamoto et al., 2013). Our simulated EMG signals are consistent with clinical EMG signals from Cernera et al. (2021), which showed EMG recordings from a patient with essential tremor during DBS-OFF and 135-Hz Vim-DBS (Cernera et al., 2021) (Supplementary Figure S6).

The EMG simulation from the model is too slow for practical implementation. Thus, we estimated the model-simulated EMG with a polynomial fit to facilitate the computational speed in a closed-loop control system. The model-simulated EMG and polynomial-estimated EMG are denoted as “reference EMG” and “estimated EMG,” respectively (Figure 5). We compared the reference EMG and estimated EMG in response to different frequencies of DBS (Figures 5, 6).

In the time domain, the estimated EMG is similar to the reference EMG across different DBS frequencies (10–200 Hz), in terms of both amplitude and variation (Figure 5A). The correlation between the reference and estimated EMGs is generally above 0.3, representing some positive correlations (Figure 5B); the correlation is not very high because of the existence of white noise in the simulations. The power is generally similar between reference and estimated EMGs (Figure 5C). In addition to the comparison in the time domain, we also compared the reference and estimated EMGs in the frequency domain (Figure 6). At each frequency in the band [2–200] Hz of the EMG activities, we computed the corresponding frequency power with Eq. 8 in the initial T = 5 s of the EMG (Figure 6). The estimated EMG [by a 25-order polynomial (Eq. 6)] is well-fitted to the reference EMG in the frequency domain, with R ² = 0.745 (Supplementary Figure S3), computed based on the signals across different DBS frequencies. Additionally, we observed other similarities between the estimated and reference EMGs in terms of the amplitude and pattern of different frequency powers of EMG activities (Figure 6). The EMG power is high with DBS-OFF and during low-frequency (<100 Hz) DBS and is mostly suppressed during ≥ 130 Hz DBS (Figure 6). During 10–80 Hz DBS, in both estimated and reference EMGs, we observed that the power is high at the harmonics of the DBS frequency (Figure 6). This might indicate that DBS could induce synchronized activities during low-frequency DBS (Florin et al., 2008; Pedrosa et al., 2013). The similarities between the reference and estimated EMGs—in both the time domain and frequency domain–indicate that the estimated EMG is a proper substitute for the reference EMG in controlling the frequencies of Vim–DBS for treating essential tremor.

We computed the system output to be implemented in a closed-loop controller updating the DBS frequency. The system output during Vim–DBS was computed with the estimated EMG (Eq. 6 ) that facilitates the implementation speed. For DBS frequency u , we defined the system output z u as the mean power spectral density (PSD) in the initial interval T = 5 s of the estimated EMG in response to DBS with stimulation frequency u (Eq. 9). PSD represents the band [2–200] Hz of EMG activities. The system output in response to different frequencies ([10–200] Hz) of DBS is presented in Figure 7:

In general, the system output decreases with increase in the DBS frequency (Figure 7, Supplementary Table S3). During low-frequency ( ≤ 50 Hz) DBS, the system output is not reduced much compared with the DBS-OFF ( u = 0) situation (Figure 7). The system output is low during high-frequency ( ≥ 100 Hz) DBS and is close to minimum during ≥ 130 Hz DBS (Figure 7). These responses of the system output to different DBS frequencies are consistent with clinical observations of the effectiveness of different frequencies of Vim–DBS (Pedrosa et al., 2013; Earhart et al., 2007; Vaillancourt et al., 2003; Dembek et al., 2020).

A PID controller (Eq. 10) was implemented to update the DBS frequency in the closed-loop system, based on the system output z (Eq. 9). The parameters K p , K i , and K d of the PID controller were chosen to be 10³ (Hz/mV²), 10⁵ Hz/(mV² ∙ min), and 5 × 10³ (Hz ∙ min/mV²), respectively; parameter tuning was performed to increase the efficacy of the controller (Supplementary Figure S2 and Supplementary Note).

As a test of the controller, when the control target is the power of the biomarker (estimated EMG, Eq. 6) from 130-Hz DBS, the PID controller can converge to the target in 10 min (Figure 8A and Supplementary Table S4). This shows that our control system is potentially effective and efficient for clinical implementation. Note that during the PID control, only the steady-state DBS frequency (reached after ∼10 min) is delivered to the patient. As we change the control target of the system output, the result of the PID control is also robustly and flexibly changed (Figure 8B). As shown in Figure 8B, the five control targets of the system output correspond to the biomarker power from both observed and unobserved DBS frequencies (Supplementary Table S4). The observed DBS frequencies (10, 50, 80, 100, 120, 130, 140, 160, and 200 Hz) were used in fitting the polynomial coefficients (Eq. 6), and the unobserved DBS frequencies are arbitrary.

The system output used in our closed-loop system is related to the power of the model-predicted EMG signals, and the optimal DBS frequency is obtained by bringing the system output to a specified control target. In clinical studies, the power of EMG is a commonly observed indicator for different movement disorders, e.g., PD (Zhang et al., 2017), ET (Halliday et al., 2000), and akinesia (Bisdorff et al., 1999). Tremor symptoms, characterized by the tremor amplitude and frequency, can be identified using the power of EMG (Hess and Pullman, 2012). Tremor amplitude is the primary indicator of the severity of tremors (Hess and Pullman, 2012). Tremor frequency can be used to partially differentiate disease types; e.g., the peak tremor frequency observed in the EMG of PD patients is often 3 ∼ 6 Hz (Zhang et al., 2017; Hess and Pullman, 2012) and in the EMG of ET patients is often 4 ∼ 8 Hz (Halliday et al., 2000; Herron et al., 2017). Thus, use of the power of EMG as a biomarker for ET in a closed-loop DBS is clinically relevant (Yamamoto et al., 2013; Herron et al., 2017). During the DBS control, the control target of the EMG power should be appropriate: a high EMG power indicates the insufficiency of tremor suppression, and a low EMG power can be related to akinesia (Bisdorff et al., 1999) and myasthenia gravis (Mills, 2005).

The ability to predict how different frequencies of DBS change neural and behavioral activities is the main advantage of our model-based closed-loop DBS. A controller (e.g., PID) can select appropriate DBS patterns that are biophysically relevant and clinically effective. Although the model parameters [for both encoding and decoding models (see Figure 2) are obtained based on sparse DBS frequencies, 5, 10, 20, 30, 50, 100, 130, and 200 Hz [human Vim data (Tian et al., 2023b) and non-human primate cortical data (Agnesi et al., 2015)], the encoding model can robustly predict the effect of an arbitrary DBS frequency in the continuous spectrum of 5–200 Hz DBS (see Figure 4 for some examples; see Table 1 in Tian et al. (2023a) for a test of robustness of the firing rate model). Model-based control of DBS was addressed in previous computational studies (Fleming et al., 2020; Grado et al., 2018; Liu et al., 2021). Despite the usability of these control systems for in silico explorations, the underlying models were not fitted to experimental data. More importantly, these models do not consider the physiological mechanisms of the DBS effects. In this work, our (encoding) model not only incorporates biophysically realistic dynamics of DBS-induced short-term synaptic dynamics but also provides an accurate fit to Vim–DBS experimental data (Tian et al., 2023b). Additionally, our data-driven decoding model provides a fit (see Figures 5, 6) to simulated EMG data in which physiological mechanisms of tremor symptoms were preserved (see Figure 2 for details of the simulation study). Our control system predicts the effect of DBS frequency on the power of EMG and delivers an optimal DBS frequency.

In recent closed-loop DBS control systems, machine learning methods have been developed to map biomarker features (input) to patients’ observed states (output) and could further deliver an appropriate DBS setting (Chandrabhatla et al., 2023; Kuo et al., 2018; Merk et al., 2022). Therefore, it is imperative to identify key biomarkers that need to be extracted from the LFP and EMG as they will serve as input features for a judiciously selected machine learning model. Castaño-Candamil et al. (2020) used a regression method to estimate tremor severity from electrocorticographic (ECoG) power in ET patients and adjusted the DBS intensity according to tremor severity. Golshan et al. (2018) used the wavelet coefficients of the STN-LFP beta frequency range as features and further developed a support vector machine (SVM) classifier for studying the behaviors of PD patients. Numerous prior studies have established high performance using SVM classifiers with input features such as phase-amplitude coupling (Chandrabhatla et al., 2023), Hjorth parameters (Oliveira et al., 2023), beta band power (Chandrabhatla et al., 2023), and burst duration (Merk et al., 2022). Using power densities within the beta (Kuo et al., 2018) and gamma (Yao et al., 2020) bands as features, hidden Markov models (Merk et al., 2022; Sun et al., 2020; Yao et al., 2020), SVM (Golshan et al., 2018), convolutional neural networks (CNNs) (Merk et al., 2022; Golshan et al., 2020; Oliveira et al., 2023), linear discriminant analysis (LDA) (Merk et al., 2022), and logistic regression (Houston et al., 2019) have been investigated. It was also recommended in a couple of studies that deep learning methods such as CNNs are worth investigating as they capture nonlinear temporal dynamics and waveform shape (Merk et al., 2022; Oliveira et al., 2023). For example, Haddock et al. (2019) developed a deep learning method that classified the behaviors of ET patients based on the PSD of ECoG and used the classification results to turn DBS ON/OFF.

A key improvement to existing machine learning methods is to incorporate physiological characterizations of the input–output mapping. In this work, we developed a physiological model to map the input (DBS frequency) to the output (model-predicted EMG). We then used a polynomial-based approximation to estimate the input–output map to facilitate the implementation speed of the control system. However, the polynomial method is prone to be less robust to unseen inputs, owing to its high-order terms (Kane et al., 2017; Beltrão et al., 1991). Thus, an important line of future work is to use state-of-the-art machine learning methods, particularly deep learning methods, to replace the simple polynomial-based input–output mapping. Consequently, it is important to understand key EMG features for muscle activation in Parkinson’s disease. The literature highlights sample kurtosis, recurrence rate, and correlation dimension as three specific EMG features that are responsive to changes in DBS treatment parameters (Rissanen et al., 2015). These features have been fed as input into the LDA, CNN, and SVM, where the SVM performed the best (Ruonala, 2022). In a few investigations, EMG features, encompassing frequency, amplitude, and regularity, were scrutinized (Khobragade et al., 2018; Wang et al., 2020). The signal mean and power of the peak frequency performed well as features when using a random forest model and a deep learning network for adaptive DBS (Khobragade et al., 2018). It is important to note that while simpler models like LDA are valued (Ozturk et al., 2020; Watts et al., 2020) for their interpretability in the context of DBS for PD, the limited availability of labeled datasets resulted in ambiguous success for complex models, particularly deep neural networks (Oliveira et al., 2023).

Our closed-loop DBS control is based on a physiological model that can generate an arbitrary amount of synthetic data, which can be implemented in fully training deep learning methods. The arbitrary amount of data in the training set will increase the accuracy and robustness of our future closed-loop DBS control based on physiological models and deep learning methods. Therefore, it is recommended that the efficacy of the SVM, logistic regression, LDA, hidden Markov model (HMM), random forests, and deep neural network models like CNNs be evaluated in greater detail using the abundance of synthetic data. Watts et al. (2020) performed a thorough retroactive study of various machine learning classifiers used to identify optimal DBS parameters for PD, and it was recommended to pursue machine learning in the context of adaptive closed-loop DBS for PD.

Synaptic depression can partially explain the therapeutic mechanisms of high-frequency DBS (e.g., Vim–DBS), which can stem from synaptic and axonal failure (Rosenbaum et al., 2014; Tian et al., 2023a). In the present work, we incorporated dynamics of DBS-induced short-term synaptic plasticity (STP)—characterized by the Tsodyks and Markram model (Tsodyks et al., 1998)—in our neural model. Such a modeling strategy of the DBS effect is consistent with previous works (Farokhniaee and McIntyre, 2019; Milosevic et al., 2021). Further details of the DBS effect can also be considered in neuronal simulations. For example, Schmidt et al. (2020) modeled the effect of DBS by generating a spherical electrical field that affects the potential of all neuronal elements (including soma and axon) within a certain distance from the DBS electrode (Schmidt et al., 2020). The electrical field induced by DBS can be non-spherical if multiple electrodes or directional leads are used (Steffen et al., 2020; Masuda et al., 2022). DBS can induce both orthodromic and antidromic activations of axons, e.g., in Vim–DBS (Grill et al., 2008) and STN–DBS (Neumann et al., 2023). In particular, during STN–DBS, the antidromic activation of the cortical circuitry is a key factor in changing neural dynamics (Neumann et al., 2023). We will investigate the DBS effect of antidromic activations in our future models and compare different models of the DBS mechanisms.

In our closed-loop system, the EMG power of a broad band is used as the system output to update the input DBS frequency. We used the band [2–200] Hz that covers the DBS frequencies [10–200] Hz, which induces DBS-evoked activities in our EMG simulations (Figures 5, 6). These DBS-evoked activities could be a mechanism of the ineffectiveness of low-frequency Vim–DBS (Ushe et al., 2004; Earhart et al., 2007; Pedrosa et al., 2013) and need to be suppressed in the closed-loop control system. Yet, EMG recordings during low-frequency DBS are very limited, and more such recordings are needed to fully investigate the underlying mechanisms. There have been closed-loop DBS systems controlling the tremor band [∼(2, 12) Hz] of EMG activities in ET patients (Cernera et al., 2021; Yamamoto et al., 2013). Cernera et al. (2021) showed that during DBS-OFF for an ET patient, most of the EMG power belongs to the band (2–12) Hz (Cernera et al., 2021). Thus, the control result may be similar when using the broad band [2–200] Hz, in which the power of [12–200] Hz is small, and the result will not be biased toward this relatively small power. In the future, we will perform the control of the tremor band ∼ (2–12) Hz of EMG activities and compare it with the current scheme.

We developed a model-based closed-loop control system that automatically updates the DBS frequency. Although the DBS frequency is a commonly tuned parameter in clinical applications (Merola et al., 2013), the tuning of another DBS parameter (e.g., pulse width and amplitude), or a combination of different DBS parameters, may also be clinically effective. Our closed-loop system adapts the DBS frequency because the underlying Vim network model was built to fit clinical data recorded during different frequencies of DBS. In the future, we will develop closed-loop systems that can adapt different DBS parameters, based on the corresponding new clinical data.

Our model-based closed-loop DBS control system is in the proof-of-concept stage for clinical implementations. The system will be implemented together with a constant monitoring of the EMG signal. In contrast to existing closed-loop DBS systems that update DBS parameters based on EMG signals solely (Yamamoto et al., 2013; Herron et al., 2017), before delivering DBS to the patient, our system predicts the effect of DBS frequency based on the underlying computational model. During implementation, the recorded EMG signal will be used to adjust our model predictions to personalize the model-based system and increase the prediction accuracy.

Our model included the clinical Vim–DBS data recorded in Vim neurons across different stimulation frequencies (10–200 Hz) (Tian et al., 2023b). However, experimental data were not sufficiently included in other components of our model. We incorporated non-human primate single-unit recordings in M1 during 130-Hz VPLo-DBS reported in Agnesi et al. (2015), but M1 activities in response to other DBS frequencies were not recorded. The EMG model simulation is consistent with some clinical observations, but was not further developed and validated by fitting clinical EMG data. In fact, in current experimental work on DBS, the EMG signal is usually recorded with DBS-OFF or high DBS frequencies (>100 Hz) (Herron et al., 2017; Vaillancourt et al., 2003; Rissanen et al., 2011), and they lack EMG recordings in response to a wide spectrum of DBS frequencies (e.g., 10–200 Hz). In the future, we plan to incorporate more experimental data into the further development of the model-based closed-loop DBS control, and these experimental data—in particular, cortical and EMG data—need to be recorded with different DBS frequencies from each individual subject.

It is worth mentioning that the synaptic connections between the Vim (VPLo in primate) are reciprocal and excitatory, though the projections that the M1 sends to the Vim are in different M1 laminae than the ones it receives from Vim (Stepniewska et al., 1994). Our model simplified this excitatory feedback relationship by considering the excitatory propagation effect as unidirectional, from the Vim to the M1. We fitted the simplified model to the recording from one M1 neuron during 130-Hz VPLo-DBS, and there are variabilities among the dynamics of different M1 neurons. In the future, we will develop a more detailed Vim-M1 model based on more M1 recordings. In this work, we modeled the spinal motoneurons as one population. Spinal motoneurons can be classified into two function groups: somatic and visceral (Fields et al., 1970). The M1–motoneuron synaptic projection can be modeled by pair-based STDP, which characterizes membrane potential dynamics using spike timings of both pre- and post-synaptic spikes (Morrison et al., 2008; Gütig et al., 2003). We will develop more detailed models of the M1–motoneuron circuits in the future. Our future closed-loop DBS systems will be constructed based on both improved models and deep learning methods.