Recordings collected from six test participants, previously reported in Tan et al. (2016) and Shah et al. (2018), were used in this study. Three participants were tested using their left and right hands, adding up to nine individuals’ datasets (Table 1). For practicality, recording across different hemispheres will be assumed to be independent, granted that factors, such as electrode type, physiology or medication may cause some dependencies. The participants were patients with Parkinson’s disease undergoing DBS surgery targeting the STN. Each gave informed consent to participate in the study, which was approved by the local ethics committee board. The participants were asked to apply a given amount of grip force for a given period following the presentation of a visual or auditory queue. Precise details on the data recording procedure can be found in Tan et al. (2013, 2016). Table 1 presents the number of trials and channels recorded. On average, a recording was 402 ± 105 s long and contained 40 ± 9 trials. Each grip was held for 4.1 ± 1.4 s. The ratio of grip against rest was 31 ± 7%, and the average holding force was 47 ± 27% of the maximum voluntary contraction. In each dataset, multiple channels of LFPs were recorded from lateral and contra-lateral electrodes (Table 1). The LFP signals were pre-processed following the same method outlined in Martineau et al. (2020). In summary, the prep-processing steps include detrending the recordings using a low-pass filter, removing main artifacts using notch filters, scaling the signals based on their standard deviation and resampling at the desired frequency. The dataset obtained was the following: an LFPs array x with time and electrode channel dimensions (for both hemispheres) such that x ∈ ℝ^N_k×N_ch, where N_k is the total number of samples and N_ch, the total number of channels. xwas re-sampled at a rate of 512 Hz, which was defined as the decoder input sampling frequency f_i. The decoder output y~ matched the detection target y at a rate of 16 Hz, defined as the decoder output sampling frequency f_o. Decoders typically operate in the range of 10–20 Hz (Shah et al., 2018; He et al., 2021). A base 2 number was chosen to facilitate downsampling, especially for implementing the discrete wavelet transform-based feature extractor. As shown in Figure 3, a sampling step-down R = f_i/f_o was performed such that y~⁢[k],y⁢[k]∈{0,1}⁢∀k∈{0, 1,…⁢NkR-1} with 0 marking a rest state and 1 an action state. A target lead of 125 ms was also introduced (Benabid et al., 2019). It is common practice to predict the target a few decoder update steps (in this case, two update steps) in advance to compensate for potential latencies introduced by processing, communication or even stimulator ramp-up (Tarvainen et al., 2004; Gruenwald et al., 2017; Benabid et al., 2019). Finally, 10% of the length of each recording was set aside for testing after tuning.

Five different frequency-domain and real-time extraction methods, including a bank of infinite impulse response (IIR) band-pass filters, a short-time Fourier transform (STFT), a recursive wavelet transform (RWT), a parameterized discrete wavelet transform (DWT), an adaptive auto-regressive moving average (ARMA) filter, based method, were implemented.

Filter banks have been proposed as a computationally efficient solution to extract standard electrophysiological sub-bands (Zhu et al., 2020). The method proposed by Gruenwald et al. (2017) was followed to construct a bank of fourth-order IIR Butterworth filters. The frequency band were defined as θ (1–8 Hz), α (8–12 Hz), β (12–32 Hz), γ₁ (32–50 Hz), γ₂ (50–100 Hz), γ₃ (100–256 Hz), leading to a total of N_b frequency bands. The power of each sub-band xp∈ℝNkR×Nc⁢h×Nb was estimated from the signal by measuring the variance of the sub-signals on a window length L and strides R. L was defined as a multiple of R such that:

where n_o is the window overlap number. A maximum L of 250 ms was allowed as in Flint et al. (2012), such that n_o ∈ ⟦0; 3⟧.

Secondly, the STFT is a reference time-frequency estimation method used in previous STN-LFP decoding studies (Bakstein et al., 2012; Hirschmann et al., 2017; Shah et al., 2018). The same epoching parameterization was used and the Kaiser window was chosen to parameterize the window shape through a single parameter β_K ∈ [0;15]. The STFT yielded a complete time-frequency transformation χ∈ℂNkR×Nc⁢h×(L2+1) of x. The same N_b sub-bands were extracted by averaging the correct frequency bins following the power operator χ[k,l,n]χ*[k,l,n] to obtain x_p. A drawback of the STFT is that the resolution of the transform is directly dependent on L. With a small number of or no overlaps, narrower lower bands are thus not properly represented (in particular α) and not included as features.

This problem is typically addressed by opting for the continuous wavelet transform (CWT) instead. In addition to better capturing non-stationary signal behaviors, the CWT allows for multi-resolution signal analysis. The Morlet wavelet in particular is predominantly used in the field and thus is also regularly employed in decoder architectures (Bashashati et al., 2016; Tan et al., 2019; Golshan et al., 2020). To reduce the computational cost and latency of the CWT, Ren and Kezunovic (2010) proposed a complex sinusoid wavelet function, ψ, shaped with an asymmetric fifth-order polynomial envelope:

Ren and Kezunovic (2010) demonstrated that the CWT with ψ can be recursively expressed. This RWT approach was implemented as a complex IIR filter bank. Each filter projected the signal on a wavelet basis. Four wavelets were used per octave between 1 and 256 Hz, or a total of 32 filters, yielding a time-frequency transformation χ ∈ ℂ^{N_k×N_ch×32}. The power array xp∈ℝNkR×Nc⁢h×Nb was obtained after a power operation and application of epoch averages on χ.

Moreover, the more common alternative to the CWT is the DWT (Bakstein et al., 2012; Mamun et al., 2015; Luo et al., 2018). The DWT is implemented using small and efficient filter bank which samples discrete points in the CWT. The main design difficulty associated with the DWT is the choice of a mother wavelet ψ (Mamun et al., 2015; Luo et al., 2018). Farina et al. (2007) proposed using the parameterized lattice realization of the quadrature mirror filter to shape the mother wavelet to neural signals which in turn helped increase decoding performance. Each node in the DWT decomposition tree was constructed using a three-tap filter. The quadrature mirror filter split a signal x into low x_g and high x_h components through two rotations θ1,θ2∈]π2;-π2[ (Vaidyanathan, 1993; Farina et al., 2007). To match f_o, N_l = log₂⁡(R) levels were required. The decomposition tree recursively filtered x_g at each level l, resulting in a set of wavelet coefficients Sc={xh1,xh2,…⁢xhl,…⁢xhNl,xgNl}. The power array xp∈ℝNkR×Nc⁢h×(Nl+1) was assembled by epoching each coefficient of S_c with a matching stride Rl=R2l, length L_l, and computing the sum of squares.

Last, a parametric modeling and spectral estimation methods have been regularly used in STN-LFP (Foffani et al., 2005; Mamun et al., 2015; Shah et al., 2018) and other neural decoders (Chisci et al., 2010). The recursive-least-square ARMA scheme proposed by Tarvainen et al. (2004) was here investigated. The algorithm extracted the ARMA coefficients x_c ∈ ℝ^{N_k×N_ch×(N_AR + N_MA)} from x. N_AR ∈ ⟦1, 12⟧,N_MA ∈ ⟦0, 6⟧ are the number of auto-regressive and moving average coefficients, respectively (Tarvainen et al., 2004). x_c was updated using a forgetting factor λ ∈]0,1[. Epoching served to average and downsample x_c before estimating of the power spectrum (Tarvainen et al., 2004). The spectrum was partitioned into N_b bands to obtain x_p.

In addition, the usage of Hjorth (HJ) parameters (activity, mobility, and complexity) (Shah et al., 2018; Zhu et al., 2020; He et al., 2021) was investigated. The signal first, x_d, and second derivatives were estimated using discrete differences. The parameters, h₁, h₂, and h₃, were then computed by taking the following the signal variance on epochs of length L and stride R:

where Var_LR is the epoch-variance operator. Hjorth parameters indirectly capture some characteristics of the signal spectrum through the time domain (Hjorth, 1970). Compared to most other time-frequency methods, HJ parameters are inexpensive to compute, making them desirable for real-time applications (Zhu et al., 2020; He et al., 2021).

A computational benchmark was conducted for every feature extraction method by measuring the time taken to extract features from a fix length (100 s) random signal and estimating the time taken to process a signal block of size R. Under real-time constraints, such a block is needed to update the decoder at a rate f_o. The computational time, thus, needs to be inferior to 1/f_o (or 62.5 ms). The SciPy implementation (Virtanen et al., 2019) of the STFT algorithm was used partly to serve as a reference, partly because it is well-optimized and uses the fast-Fourier transform algorithm. Given the relatively large size of the extraction windows, a direct filter implementation would have been inefficient. Aside from the STFT, all extractors were implemented using the NumPy library (Harris et al., 2020). In every extractor, every filter was recursively updated at each time step to imitate real-time conditions. When implementing the step update, the well-optimized multi-dimensional array operation capabilities of NumPy (Harris et al., 2020) were carefully considered.

From a real-time deployment perspective, it could be costly to carry signal acquisition on all channels, extract all features and only select a sub-set for decoding. Instead, a channel selection approach was chosen whereby channels can be de-activated if non-selected. Channels were selected based on the cumulated mutual information (MI) of their extracted features in relation to the decoded target. Given x_p and y, the average information content M⁢I¯ of the l^th channel was defined as:

where MI is the mutual information operator (Pedregosa et al., 2011; Ross, 2014), N_i is the total number of features for each channel, so either N_b, N_l + 1 (if DWT is used) or 3 (if HJ is used), and l is the channel index. M⁢I¯ was sorted in descending order. N_c ∈ ⟦1,N_ch⟧ top ranked channels were selected. Features of the selected channels were further processed by applying a logarithm transform, normalization, and further smoothing using the feature-wise Kalman filter proposed by Gruenwald et al. (2017) (except for the ARMA extractor, where it could be adjusted through λ). The feature matrix X∈ℝNkR×Nf, with Nf=Nc⁢Nw⁢Ni, was assembled by flattening all the channels and bands (or coefficients in the DWT case, or parameters in the HJ case) dimensions and concatenating N_w ∈ ⟦0; 15⟧ past time points to form a moving window of features (or a maximum length of 1 s) (Flint et al., 2012; Shah et al., 2018; Tan et al., 2019; He et al., 2021). The feature and label pair (X,y) were under-sampled according to repeated-edited-nearest-neighbors method to create a more balanced training set before classifier fitting (Tomek, 1976; Lemaître et al., 2017).

Furthermore, the parameters n_o, N_c and N_w controlled the total decoder window τ_d (Figure 3) and the total number of features N_f. These two metrics were used to evaluate the complexity of the resulting decoder. A decoder with a smaller τ_d established more temporally located relationships between neural features and the decoded state. Practically, decoder with longer τ_d required more memory to store past signal epochs for extraction and past feature epochs for state estimation. Likewise, N_f also directly influenced memory storage and, more importantly, the computational cost of updating the state estimator. Tuning N_c and N_w also acted as a form of feature selection, which could help reduce model complexity.

Four different classifiers retained from He et al. (2021) were tested to estimate y from X: (1) the Gaussian Naïve Bayes (NB) classifier was chosen for its lack of tuning parameters (Mamun et al., 2015). The NB served as a special case to evaluated if decoder parameters (e.g., n_o, N_c, or N_w) other than fitting regularization constants could be leveraged to increase performance. (2) The linear-discriminant analysis (LDA) classifier has been shown to deliver good decoding performance (Flint et al., 2012; Gruenwald et al., 2019). A shrinkage parameter γ_c ∈ [0,1] was used to regularize the covariance estimation of X. (3) Logistic regression (LR) has successfully been applied to STN-LFPs decoding (Shah et al., 2016; Tan et al., 2019; He et al., 2021), particularly for sparse model fitting usingL₁-norm regularization. The regularization was set through a parameter C ∈]0,5]. Its setting acted in part as a form of feature selector. (4) Thanks to their versatile non-linear kernels, support vector machine (SVM) classifiers have been regularly employed for high-performance decoding (Mamun et al., 2015; Golshan et al., 2016, 2018; Nandy et al., 2019; He et al., 2021). The radial-basis variant was chosen. It required a regulation and kernel parameter C,γ_c ∈]0,5] (Nandy et al., 2019). Given the expandable nature of the feature space (of size N_F), the computational demand of the SVM fitting algorithm and the fact that decoder fitting occurred across full-length time-series at each BO iteration, principal component analysis (PCA) dimensionality reduction had to be introduced to compress X by a factor r_c ∈]0,1] (Golshan et al., 2018). The maximum number of selected components was as well capped at 64. Once deployed, the computational cost of operating the NB, LDA, and LR scales with N_f. For the SVM, it also scales with the number of support vectors used by the classifier. Other non-linear methods such as simple decision trees and nearest-neighbor algorithm were not considered as they were found to be less performant (He et al., 2021). In all, a range of incrementally more sophisticated, and arguably more complicated to tune, classifiers was chosen.

Moreover, to denoise the output prediction of the classifier, a pair of thresholds operated as a transition trigger. Given the classifier probability p(y = 1)[k], the predicted state y~⁢[k] was given by:

where T ∈ [0,1] is the upper threshold, and r_T ∈ [0,1] is the position of the lower threshold relative to T. Thresholds are typically set at the optimal point of the receiver operating characteristic (ROC) curves, balancing the sensitivity and specificity of a classifier. Unfortunately, this curve cannot be drawn for a dual-threshold system. Instead, the geometric mean g between sensitivity and specificity was chosen to estimate an optimal trade-off (Cai et al., 2010). Given the detection true-positive rate (TPR) and false-positive rate (FPR), g was given by:

If TPR = 0, then g = 0. Likewise, if FPR = 1, then g = 0. Consequently, g only leans toward 1 if both objectives are balanced, thus avoiding the extreme situation where, for instance, the state-estimator only detects true positives or negatives. On a traditional ROC diagram, g would indicate a curve point sticking furthest toward TPR = 1 and FPR = 0. An eightfold cross-validation was employed to evaluate the average detection performance g¯. The number of folds was chosen to be double that of the number of cores available on the processor (Intel^® Core i7-7700K 4.20 GHz) to accelerate the evaluation of the folds through parallelization.

Figure 1 illustrates different parameters introduced in the different processing stages and, Figure 2 the BO procedure. BO was implemented to optimize g¯ using the scikit-optimize library (Head et al., 2021). The parameters formed a bounded search space S from which a parameter vector ρ ∈ S was sampled. To quantify the performance increase resulting from tuning, initial quantities, reflecting values typically found in other studies, formed a default parameter vector ρ₀. A 3-tap Daubechies wavelet was used for the DWT by default (Rieder et al., 1998). A summary of all tuned parameter, search intervals and default values is given in Table 2. Given g¯ as a function of x, y, and ρ, the best parameters ρ*were defined as:

The corresponding maximum cross-validation performance was noted as g¯*⁢(ρ*,x,y). The BO scheme repeatedly fitted and evaluated the pipeline at each iteration (until 90 iterations were reached). The algorithm was seeded by first evaluating ρ_o. In a “tell” step, the optimizer fitted a surrogate model which mapped S to g¯. In an “ask” step, the optimizer, following the policy of a given acquisition function, suggested a new ρ (Shahriari et al., 2016). A Gaussian-process with Matern Kernel (Bashashati et al., 2016) was chosen for the surrogate model alongside the Gaussian-process hedge acquisition function (Shahriari et al., 2016; Head et al., 2021).

Before optimizing every feature extractor and classifier pair, a preliminary comparison between BO and RD was conducted by optimizing the STFT-LR decoder for three participant recordings (PD01R, PD05R, and PD06L), which represent well different optimized performance levels within the population sample. The optimization procedure was, in this case, repeated five times with different random seeds to produce average convergence trends. Individual trends were obtained by noting for every i-th iteration of the maximum performance obtained on or before said iteration g¯i*. Average cumulated regret trends were also produced for both algorithms (Shahriari et al., 2016). A regret trend r_i was defined as:

where g¯i is the score obtained on the i-th iteration and N_i, the total number of iterations. An initial 60 iterations were set. g¯i*-g¯i measures the quality of every decision the optimizer takes. Perfect decision-making would find a new optimum point at every step such that g¯i*-g¯i=0. Both convergence and regret trends gave insight into the progression of the optimization sequence.

Following the preliminary comparison, the gathered evaluation criteria summarized in Figure 2 were collected for every feature-extractor-classifier pair. The maximum cross-validation performance g¯*, TPR and FPR and the optimal parameter vector ρ* were recorded. The significance of the absolute improvement g¯*-g¯⁢(ρ0,x,y) was evaluated using a Wilcoxon paired test. A Mood test was also performed between g¯* and g¯⁢(ρ0,x,y) to measure any reduction in inter-subject variance. Performance evaluation also included plotting the decoder output prediction on a test segment. Sensitivity analysis of the BO models was used to relate the contribution of each parameter to g¯* using the SALib library (Herman and Usher, 2017). S was sampled using the Saltelli method (Saltelli, 2002), and the BO surrogate model was analyzed using the Sobol method (Sobol, 2001). The Sobol indexes measured the ratio of individual parameter input variances to surrogate model output variance (Sobol, 2001). Total order indexes for every parameter were recorded to capture both first-order and higher-order interactions (Herman and Usher, 2017). Analysis of variance (ANOVA) and post-hoc Tukey pairwise tests were applied to measure the dependence of g¯⁢(ρ*,x,y), TPR, FPR, improvement, ρ*, τ_d, N_f and parameter sensitivity indexes on the choice of extractor and classifier. Statistical testing was conducted using the SciPy (Virtanen et al., 2019) and Statsmodel (Seabold and Perktold, 2010) toolboxes. Figures were produced using the Seaborn (Waskom, 2021) and Matplotlib (Hunter, 2007) toolboxes.

Figure 4 presents averaged convergence and cumulated regret trends computed for three representative participant recordings. The trends indicated that the RD scheme was very effective at uncovering good initial parameter vectors, whilst BO had a rather slower convergence rate starting from ρ_o. Nonetheless, after a few iterations, the convergence rate and plateau of the BO algorithm overtook that of the RD scheme. The cumulated regret trends confirmed that BO mapped sufficiently well the parameter space with each iteration, thus decreasing the g¯i*-g¯i through informed decision-making. As expected, since RD is memoryless, that gap remained constant, leading to a linear cumulation of decision errors.

Thus, for the presented decoder tuning problem, BO was better at ordering good solutions (i.e., g¯i is more likely to be superior to g¯i-1, which is important for sequential learning) and at fine-tuning in later iterations. Consequently, for the remainder of the tuning runs, to bootstrap the BO search starting at ρ_o, five additional random samples were drawn to seed the initial fit of the BO surrogate model. Thirty more iterations (adding to a total of N_i = 90 iterations) were added to allow for further fine-tuning.

Figure 5 reports the optimal cross-fold average g¯⁢(ρ*,x,y), TPR and FPR achieved after BO tuning for every extractor and classifier pair. Medians laid within 0.6–0.75 for the g¯ and TPR, and within 0.20–0.35 for the FPR. Combinations of ARMA STFT extractors with LDA and LR classifiers achieved the highest average g¯ performance of 0.74 ± 0.06 (mean ± SD). ANOVA testing on g¯ suggested that the choice of classifier and extractor impacted performance (p < 0.01). Post-hoc testing showed that all classifiers and extractors performed comparatively well except for the NB (p < 0.05) and the HJ (p < 0.05). The TPR depended only on the selected extractor (ANOVA, p < 0.01), with the HJ underperforming relative to the ARMA, IIR, and STFT (post-hoc tests, p < 0.05). Conversely, the FPR depended on the choice of classifier (ANOVA, p < 0.05), with the NB underperforming relative to the LR (post-hoc tests, p < 0.05).

Figure 6 reports the improvement in the overall decoding performance with the BO-based parameter optimization compared to that with default often used parameters fixed across all participants [g¯⁢(ρ*,x,y)-g¯⁢(ρ0,x,y)]. Paired testing showed that this gap was significant in every case (Wilcoxon, p < 0.05). ANOVA testing showed that the improvement gap depended again on the classifier and extractor pair (p < 0.01). Post-hoc testing indicated that improvement margins were tighter for the LR and LDA classifiers (p < 0.05). Default parameter solutions for the SVM-PCA and NB led to poorer decoding solutions. Lastly, no significant reduction in inter-subject variance was measured between default and optimized parameter performance (Mood, p > 0.05).

Figure 7 presents the best-decoded test performance for every user. The classifier output signal p[k] varied noticeably in terms of quality and temporal dynamics throughout the cohort. Effective detection occurred with high and sufficiently sustained p throughout the action with manageable levels of baseline noise (when y[k]=0) and action noise (when y[k]=1), such as observed in PD03L, PD04L, PD02L, and PD01R. In other instances, p adopted a more impulsive behavior, such as in PD02R or PD02L, making a full capture of the action difficult. Although the upper threshold prevented high baseline noise from triggering the detection, such as in PD03L, PD04L, or PD01R, it could not contend with very large spurious impulses outside of the action (when y[k]=0) as those found in PD04L, PD05R, PD06R, and PD06L. Conversely, spurious drops would also occur within the action (when y[k]=1) such as in PD05L, triggering the lower threshold.

Optimal parameter distributions were tested using ANOVA. The test results suggested that the thresholds (T and r_T) did not depend on either the classifier or the extractor (p > 0.01) but rather on the intrinsic signal properties of individual recordings. The number of channels (N_c) depended on the choice of classifier (ANOVA, p < 0.01) but not of the extractor. Post-hoc testing suggested that fewer channels were selected using the NB classifier rather than the LR or LDA (p < 0.05). The number of overlaps between windows (n_o) was dependent on both estimator and extractor choices (ANOVA, p < 0.01). The STFT adopted a significantly larger extraction window length (L) than the ARMA, IIR, and HJ extractors; the LR, than the LDA and NB classifiers (post-hoc tests, p < 0.05). Likewise, the choice of feature window length (N_w) depended on the extractor (ANOVA, p < 0.01) but not on the classifier selection. Notably, the STFT, IIR, and HJ required significantly larger N_w than the ARMA extractor (post-hoc tests, p < 0.05). Figure 8 presents distributions for the decoder analysis window (τ_d) and the total number of features employed for decoding (N_F, section 2.3). τ_d and N_f depended on the choice of extractor and classifier (ANOVA, p < 0.01). Because of the smaller optimum N_w and n_o, the ARMA extractor and the NB classifier required significantly smaller windows than the STFT extractor and LR classifier (post-hoc tests, p < 0.05). Likewise, because of smaller N_w, N_c or N_i, the decoders using the ARMA and HJ extractors or the NB classifier required significantly fewer features than the STFT extractor or LR classifier (post-hoc tests, p < 0.05). Lastly, there was no dependence between optimal extractor parameters and classifier choice and no dependence between classifier parameters and extractor choice.

The index distributions yielded from the sensitivity analysis are presented in Figure 9 for each parameter. Considering the subset of common parameters (n_o,N_c,N_w,r_T,T), ANOVA testing suggested that the importance indexes varied depending on the parameter type and the classifier selected (p < 0.01), but not on the extractor (p > 0.05). Post-hoc testing showed that the upper and lower thresholds in the final decoding stage (T and r_T) were significantly more important than all other parameters (p < 0.05). The window overlap (n_o) was significantly more important than the number of channels selected (N_c) and the window length (N_w) (p < 0.05), despite clear preferential settings depending on the selected classifier or extractor. Post-hoc tests also showed that the NB and the LDA attributed more importance to those parameters than the LR classifier (p < 0.05). Post-hoc tests were performed for extractor combinations with the LR classifier with all available parameters. Aside from T,r_T, the regularization constant, C, was the most important parameter (p < 0.05), with the unique exception of the ARMA λ parameter. In contrast, with the NB classifier, all parameters had shared similar levels of importance (p > 0.05), except for n_o, which was more important with ARMA and IIR pairing (p < 0.05). The remaining parameters had comparable importance with the LDA and SVM-PCA classifiers. Regularization parameters did not hold as much importance as they did with LR classifiers. Lastly, extractor-specific parameters other than λ (the model order, N_AR,N_MA, for ARMA, the lattice parameters, θ₁,θ₂ for DWT and the window shape factor, β_K, for the STFT) were never significantly more important than other parameters (omitted some exceptions, Figure 9).

Figure 10 presents a benchmark to evaluate the computational time taken to process features for a single decoder step for each method. All of them were sufficiently fast to operate in real-time. A large overhead was observed when features were extracted from a single channel, most likely associated with slow array manipulation in python at each time-step. More complex methods such as ARMA and DWT were more penalized than others (more linear transformation for the ARMA and multi-level processing for the DWT extractor). The DWT should be, in practice, far more efficient (Farina et al., 2007). Nonetheless, extraction processing time was found to scale favorably as multi-dimensional signals were introduced. More demanding processes, such as the RWT, ARMA and, to an extent, the IIR extractors, had a notable increase in computational time with every added channel. Other methods, such as the DWT and HJ were more computationally efficient. The benchmark proposed in this study is specific to the programming language, operating system and processor. Different hardware and software optimization considerations should ultimately be made based on the target deployment device and the selected method.

The optimum performance obtained in this study is comparable to the results of past studies that have employed similar asynchronous and real-time constraints to decode STN-LFPs (Mace et al., 2014; Niketeghad et al., 2018). The limitation remains that stable detection can vary significantly across participants. We note, however, that BO did not reduce inter-subject variance, indicating that lower-performing cases may be due to a lack of information related to the decoding target in the recorded signals. Moreover, the performance increase can be achieved from the simultaneous optimization of the extractor and classifier, indicating that default parameters may not capture a priori expertise knowledge. A clear benefit of this approach was observed principally in the NB classifier. Nonetheless, end-to-end decoder optimization may not have been the most efficient approach in the proposed architecture and decoding task context. The tuning problem could potentially be broken down into several stages, e.g., one for feature extraction (optimizing some desirable attributes of X) and another for state estimation (optimizing g¯). Separating these two stages would also be more practical to implement estimators with more computational intensive fitting such as the kernel SVM (and avoid resorting to PCA to compress X).

Furthermore, the decoding of movement in this investigation provides a basis to apply the method to other DBS targets, state types (e.g., tremor, gait, etc.), and BCI tasks to further gauge the impact of tuning. Online testing is needed to perform decoding sequentially with stimulation in an adaptive configuration (He et al., 2021), continuous tuning (Grado et al., 2018), and for deployment, accounting for the fixed computational and memory budget of an embedded device (Zhu et al., 2020). Finally, findings in this study can further guide algorithm selection, and a BO approach could help ease deployment by rethinking the task objective with these budgeting constraints in mind.

The optimization process involved a range of parameterization trade-offs. For instance, channel ranking or other decomposition strategies are not a strict guarantee of improved performance through a reduction in classifier overfitting (Shah et al., 2018). Except for the NB, the BO algorithm kept adding channels despite diminishing returns in performance gain (Flint et al., 2012). A multi-objective approach or optimization constraints could narrow the set of reasonable tuning solutions. Penalty for delayed detection, output chattering, or model complexity could be introduced (Mace et al., 2014). Parameters such as the number of channels selected (N_c) and the window length (N_w) could then become more relevant. Further practical improvements could also be made to the BO algorithm implementation by introducing robust convergence criteria (Bashashati et al., 2016) and more adaptive acquisition functions (Grado et al., 2018).

For feature extraction, parameters influencing the temporal smoothness of the feature signals, such as the number of overlaps between feature extraction windows n_o and the ARMA forgetting factor λ were important to decoding. A second-order Kalman filter (Chisci et al., 2010; Yao et al., 2020) could be tuned using BO and applied either in the feature space for further smoothing or at the classifier output to help stabilize the state transition. Reducing the decoder output frequency f_o could also help reduce higher-frequency noise at the cost of a slower refreshing rate. An intermediate cost-function measuring signal-to-noise ratio could be employed to tune extractor parameters in the first optimization stage. Feature extractors overall delivered comparable levels of performance. With hindsight, testing other signal features that were not purely spectral-based (Bakstein et al., 2012; Wang et al., 2018) or different combinations of feature bands might have been more intriguing (Tan et al., 2016). The BO algorithm could have been used to select the combination of features better suited to decoding.

Although the extractor was not critically associated with the detection performance, this work suggests that the choice of extractor still influences extraction and feature window lengths and the number of features used for classification, elements important to the decoder design. Some extractors may thus be better suited to some decoding tasks than others. For instance, the ARMA extractor was among the more performant and instantaneous methods. It can also be helpful if connectivity features (Mamun et al., 2015) or instantaneous phase estimation (Foffani et al., 2005) are of interest. However, it comes at the cost of a specialized implementation and a more computationally intensive process. The STFT, in comparison, is a practical, efficient and versatile extractor. Since the window shape only had a limited influence on decoding, using, for instance, a recursive formulation of the exponential window (Grado et al., 2017) could help minimize windowing lag. Although, the performance with the RWT did not supersede that of the STFT, it has the benefit of combining a narrow-band filter in a compact IIR form with a complex phase estimator, which could be useful to isolate a single or few specific sub-band phase and power. Otherwise, the additional computational burden of a large number of relatively long complex filters associated with a full-signal decomposition does not seem justify. Smaller filter banks, such as the IIR extractor, are relevant to embedded applications (Zhu et al., 2020). The mutual information-based channel selection scheme could be extended to incorporate only informative power bands. Selected individual features across contacts is here more practical as only the relevant filters can be activated on a given channel greatly reducing the computational burden. Given the comparatively short τ_d for the DWT and its computational efficiency, BO could permit the optimization of longer wavelets which could help match the performance of other extractors. The packet transform could introduce a sparser and more adaptive signal representation (Mamun et al., 2015; Luo et al., 2018). Lastly, the HJ parameters were not sufficiently informative and should only be used as complementary features for movement decoding (Shah et al., 2018).

The double-thresholds system played an important role in providing stable state transitions. It was observed that the upper threshold was generally set above spurious events to prevent false-positive triggering. A higher triggering level delayed detection, especially if the transient response of the feature extractor and classifier combination was slow. A lower second threshold, close to the baseline activity of the classifier output, was preferable. It aligned the detection and the action ending better and prevented a spurious transition if the classifier output suddenly dipped during an action. BO could thus help calibrate more complex state machines required for adaptive DBS (He et al., 2021).

Aside from the NB classifier, there was no substantial difference between the classifier types in terms of performance. However, not every classifier was best suited for the proposed detection task. Fitting the SVM-PCA classifier was more computationally demanding, given the temporal longitudinal aspect of the dataset. Adjusting f_o could help in this regard (Golshan et al., 2016, 2018, 2020). In addition, the SVM algorithm does not yield prediction probabilities. Given the importance of thresholding for state transition, those cannot be omitted and thus have to be estimated in an additional post hoc process (Platt, 1999). LR and the LDA classifiers may be thus preferable with the LR was more easily calibrated through regularization. Ensemble methods, using a set of classifiers rather than a single unit (in particular decision trees), have proven to outperform other machine-learning methods when deployed on small size, complex, and noisy BCI datasets (Gruenwald et al., 2019; Zhu et al., 2020) and could be easily implemented using the proposed method. The method could also be applicable to neural-network based decoders (Ahmadi et al., 2019; Golshan et al., 2020) specifically tuning of decoder specific parameters such as selected channels, thresholds or window convolution size.