Our main interest was the initialisation, development and the average force during the ‘holding phase’ of each grip; therefore we focused on the period of time from one second before the cue to 2.8 s after the cue (before force releasing) for decoding. The force trajectory of each individual trial of each subject was normalized against the average maximal force that subject achieved in their maximal effort trials. We hypothesised that the relationship between LFP features and generated grip force (the transfer function) could be captured by a first order linear dynamic model (Equation 1):

Where ${\displaystyle K_p}$ is the steady state proportional gain, $T_p$ is the time constant of the first order system which is a measure of how fast the force output responds to brain signal changes, and $T_d$ is the time delay between the brain control signal (LFP) and the measured force which describes the latency between the timing of brain signal changes and force output onset. Different models with different assumptions about what is the effective force control feature from the STNr (STN region) LFP and how different features are combined to encode force were tested (see Table 2 and details in Materials and methods).

Cross-validation, within each STN/contralateral hand as well as across different patient groups, was used to evaluate the generalisability and decoding performance of the proposed models for predicting gripping force (see details in Materials and methods). Within each STN/contralateral hand, after the model parameters had been fitted based on training data, several key variables were quantified to evaluate the performance of the model in predicting force in another set of test data recorded from the same hand: (1) the correlation coefficient between the predicted force and measured force (WithinTrialR), which can be used to evaluate the prediction accuracy of the force development and force trajectory within each individual trial. (2) The root mean square error, which quantifies the average distance between the prediction and actual measurements normalized by the maximal measured value (nRMSE). (3) The correlation coefficient between the average predicted force at the holding phase (over 1–2 s after cue onset, when the gripping force was relatively stable) and the measured force at the holding phase across different trials (StableFrcR). This was used to evaluate the prediction accuracy of the stable target force across different effort levels. (4) The difference between the timing of the predicted force onset and the timing of the measured force onset (DifRT).

The performances of different models in predicting force averaged across multiple trials using different STN features (see Materials and methods for details) were first evaluated. To achieve this, gripping force and LFPs measured from the contralateral STNr were grouped into low effort trials (with self-reported effort <= 5, trial number = 15±1) and high effort trials (with self-reported effort > 5, trial number = 14±1). Average STNr LFP features and force trajectories were calculated for each effort condition. The STNr LFP features and force for one effort condition were used to estimate the model parameters (model fitting), and the models were then used to predict force for the other effort condition (model testing). The within-trial correlation coefficients between the predicted force and measured force (WithinR) and the nRMSE of the predicted force from different models during model testing are shown in Figure 2A and B. There was no significant difference consistent across the STNs between the three models using beta and gamma ERS as model inputs (Models 1–3) in the predictive performance in terms of either WithinR or nRMSE. These three models with both beta and gamma activity as inputs performed better than models in which extra information about alpha activity was also included (Models 4,5), or models with activities from a single frequency band as model input only (Models 6–8). Figure 2C shows the BIC values combining the force prediction for low effort and high effort for all models. One-way repeated ANOVA identified a significant effect of models in the BIC (F(7,56)=11.777, p<0.001). Paired t-tests showed that Model 2 was significantly better, in terms of BIC, than Models 4–8 after FDR multiple comparison correction: (t(8)=−3.63, p=0.007 compared with Model 4; t(8)=−4.950, p=0.001 compared with Model 5; t(8)=−9.347, p<0.001 compared with Model 6; t(8)=−5.757, p<0.001 compared with Model 7; t(8)=−3.675, p=0.007 compared with Model 8).

Figure 2

Download asset Open asset

Figure 2—source data 1

The excel data file related to Figures 2 and 6.

It contains data about median WithinR and the normalized RMSE when different models were used for force prediction, as well as the average power modulations at the beta band and gamma band for each STN that was recorded in the main paradigm and used for the analysis.

https://doi.org/10.7554/eLife.19089.008

Download elife-19089-fig2-data1-v2.xlsx

Figure 3 shows how Model 2 with beta and gamma having different linear gains can be used to predict the force generated by the contralateral hand for individual STNs. In Figure 3A, data averaged across low effort trials (with self-reported effort <= 5) were used for the estimation of model parameters (model fitting). The predictive performance of the model was evaluated on the data averaged across high effort trials for the same STN (with self-reported effort > 5, Figure 3B). Figure 3C shows the fitting and predictive performance of the model on data from the other 8 STNr electrodes in which consistent movement-related power modulations in both the beta and the gamma band were observed. The correlation coefficient between the predicted force and the measured force across all STNs was between 0.902 and 0.987 with a median value of 0.952 for high effort force (data during low effort used for model parameter estimation) and between 0.921 and 0.992 with a median value of 0.969 for low effort force (data during high effort used for model parameter estimation). The nRMSE of the prediction force was between 5.60% and 28.07% with median value of 16.5% for high effort levels and between 10.5% and 23.8% with median value of 16.4% for low effort levels. The differences in reaction times (DifRT) between the predicted force and the measured force were 11 ± 27 ms for high effort levels and −26 ± 31 ms for low effort levels. These were not significantly different from zero (t = 1.25, p=0.247 for high effort levels and t = 1.03, p=0.333 for low effort levels, one-sampled t-test compared to zero). The estimated Td ranged between 0 and 264 ms for low effort trials and ranged from 0 and 277 ms for high effort trials across STNs. Together, the results demonstrated that a simple first order linear dynamic model with beta and gamma as inputs with different linear gains can be used to describe the relationship between STNr activity change and measured force; and that, in addition, the model with parameters derived from one set of data can be used to predict force exerted at other effort levels in the same subject.

Figure 3

Download asset Open asset

Figure 3—source data 1

The Matlab data file containing source data related to Figure 3.

The variable names are self-explanatory. ‘DataPlotC.Axisi’ contain data for the ith axis in PlotC. Column 1–7 in this variable corresponds to time, predicted force for a low effort, predicted force for high effort, fitted the force for a high effort, fitted the force for a low effort, measured the force for a high effort and measured the force for a low effort, respectively.

https://doi.org/10.7554/eLife.19089.010

Download elife-19089-fig3-data1-v2.mat

Figure 4 shows the performance of the model with beta and gamma power change in the STNr LFP as inputs (Model 2) in predicting contralateral gripping force in individual trials for one exemplar STN. For half of all the individual trials (n = 21) recorded from this STN and contralateral hand, the within-trial correlation coefficient between the predicted force and measured force was equal or larger than 0.78. The actual measured stable force during the holding phase and the predicted stable force were then quantified for each individual trial (Figure 4D). Figure 4D shows that the stable force during the holding phase varied from zero to 100% of the maximal voluntary force across all the trials. Linear correlation was applied to the actual measured stable force during the holding phase and the predicted stable force of all individual trials. This predicted stable force correlated with the measured stable force (n = 21, r = 0.815, p<0.001), suggesting that the STNr LFP in conjunction with the dynamic model can predict both the force trajectory in individual trials and the stable force achieved across trials with different effort. The line of best fit between the predicted stable force and measured stable force was y = 0.96x + 2.4. The regression gradient close to unity suggests that the prediction matches well with the measurement with no systematic overestimation or underestimation. During one trial for this patient, there is a force increase predicted from the STNr LFP, but there is no actual measured force change in the dynamometer (indicated by ** in Figure 4B). This may be due to changes in the STNr LFP signal following the cue without the actual force change being registered in the dynamometer, as might arise when movements of the limb did not involve the dynamometer.

Figure 4

Download asset Open asset

Figure 4—source data 1

The Matlab data file containing source data related to Figure 4.

The variable ‘DataPlotA’ is a data structure containing fields about the movement-related power changes over time (‘ERS’) for one exemplar subject. ‘DataPlotB’ is a data structure containing fields about measured force (‘MeasForce’) and predicted force based on Model 2 (‘EstForce’). ‘DataPlotD’ is a data structure containing measured stable force and predicted stable force for each individual trial and the linear regression between them.

https://doi.org/10.7554/eLife.19089.012

Download elife-19089-fig4-data1-v2.mat

Considering all the STNs in which significant movement-related modulations were observed in both beta and gamma bands (N = 9), the correlation between the predicted stable force and measured stable force across trials ranged from 0.448 to 0.913 for the contralateral hand when Model 2 was used, with the regression gradient between the predicted stable force and measured stable force ranging from 0.91 to 1.12 (not significantly different from unity across STNs based on one-sampled T-test, p=0.431). There was no significant difference between the three models which took into account both beta and gamma activities (F(2,16)=0.207, p=0.815, Figure 5A). In addition, the distribution of the within-trial correlation coefficients across all trials was not significantly different between Model 2 and Model 1 (ks-test p=0.2850), Model 1 and Model 3 (ks-test p=0.9584) or Model 2 and Model 3 (ks-test p=0.9228). For convenience therefore, Model 2 was taken as representative of Models 1–3. Considering all the individual trials from the 9 STNr, the median value for the within-trial correlation coefficient was 0.732 for Model 2 (Figure 5B and C). There were around 13% of trials during which the correlation between the measured force and the predicted force were negative, suggesting a failure in predicting force trajectory.

Figure 5

Download asset Open asset

Figure 5—source data 1

The Matlab data file containing source data related to Figure 5.

The variable ‘WithinR_Test_Contr’ contains data for the WithinR for each individual trial when Models 1–5 were used to predict force generated by the contralateral hand, for the 9 STNs in the main paradigm that are included in the main analysis. ‘WithinR_Test_Ipsi’ contains the WithinR when the models were used to predict force generated by the ipsilateral hand. ‘StbR_Test_Contr’ and ‘StbR_Test_Ipsi’ contain the correlation coefficients between the predicted stable force and measured stable force for the 9 STNs when different models were used.

https://doi.org/10.7554/eLife.19089.014

Download elife-19089-fig5-data1-v2.mat

The predictive performance of the STNr LFP for the force generated by the ipsilateral hand was also evaluated. The cross-trial correlations between the predicted stable force and measured stable were significantly lower for the ipsilateral hand than for the contralateral hand (F(2,16)=10.629, p=0.012). The distributions of WithinR for ipsilateral force prediction were significantly different from that for the contralateral force prediction (ks-test p<0.001 no matter whether Model 1, Model 2 or Model 3 was used), indicating that the STNr LFP was more informative in predicting force generated by the contralateral than the ipsilateral hand. However, the degree of lateralisation could have been underestimated due to the potential presence of mirror movements, even though we asked patients to avoid these and in half those recordings used for decoding we also had bilateral upper limb EMG recordings.

Figure 6 shows how the performance of force prediction based on the STNr LFP changes with movement related reactivity in the gamma band and beta band. Here we considered all the STNr with significant movement-related modulations, whether either in one or other, or both frequency bands of interest (N = 12). The median values of the within-trial correlation coefficients and stable force correlation coefficients increased with the average movement related synchronisation in the gamma band, with an exponential model $y=a*e^{-bx}+k$ explaining 75.7% (p<0.001) of the variance in WithinR and 88.1% (p<0.001) of the variance in StableR across STNs. There was also a trend for better prediction of force with increasing movement related beta desynchronization, but linear regression fitting was not significant for either WithinR or StableR. If we assume that movement-related synchronisation in gamma activity upon gripping is a good proxy for proximity to the motor region of the STN (and perhaps the upper limb representation within it) then these findings suggest that the electrode has to be very near to this region if recorded activity is to have a decent prospect of force prediction. This assumption was borne out by the fact that the movement-related synchronisation in gamma activity dropped by 71 ± 7.7% from the bipolar channel showing the most movement-related modulation to the average modulation in the remaining two bipolar channels. This drop was not so acute, 59 ± 7.9%, for movement related modulation in the beta band, which also only showed a trend for better prediction of force as movement related modulation in this band increased. Finally, there was a lack of significant correlation between the accuracy of the prediction of force and baseline beta (Spearman r = 0.468, p=0.13 for WithinR; Spearman r = 0.363, p=0.25 for StableR) or gamma activity (Spearman r = 0.281, p=0.37 for WithinR; Spearman r = 0.273, p=0.39 for StableR) measured during rest across subjects, suggesting that it is the reactivity of the power changes during gripping that is important for prediction, perhaps because it is more specific for the corresponding motor representation.

Figure 6

Download asset Open asset

Finally, the first order dynamic model with beta and gamma power changes as inputs (Model 2) was used to predict the gripping force from individual trials in an independent patient group that performed a different, but related, paradigm. Ten patients were recorded in this study and were asked to grip as fast and as strongly as they could in each trial in response to an external cue. Twenty trials were collected for each hand when the patients were on their normal dopaminergic medication. Details of patient information and experimental paradigm were previously reported (Anzak et al., 2012). The average movement related power change in the STN LFP activities and results of force prediction in this patient group are shown in Figure 7. The median correlation coefficient between predicted force and measured force for individual trials (median WithinR) ranged between 0.3848 and 0.9421 for the 20 STNs and contralateral hand (Figure 7B). Considering all the 397 individual trials across all STNs, more than 50% of the trials had WithinR more than 0.7859 (Figure 7C and E). We also found that incorporating alpha activity (Model 4) improved the force prediction accuracy for maximal effort gripping in this patient group. When Model 4 was used, the median WithinR ranged between 0.6731 and 0.9625 across all STNs, and more than 50% of all individual trials had WithinR more than 0.8699. BIC analysis showed that Model 4 (BIC = 2116.2±135.7) was significantly better than Model 2 (BIC = 1873.0±110.3) for this patient group performing maximal effort gripping (∆BIC = 243.3±95.9, t(19)=2.5369, p=0.020 comparing the BIC values, Figure 7D).

Figure 7

Download asset Open asset

Figure 7—source data 1

The Matlab data file containing source data related to Figure 7 (data from an independent patient group on a maximal effort gripping paradigm).

The variable ‘ERS_AllSTN’ contains the average movement related power changes from 1 to 90 Hz in the LFP signal from each STN recorded. ‘WithinR_IndividualSTN’ contains data of WithinR for each individual trial recorded from each individual STN when Models 1–5 that are tested on this paradigm. ‘WithinR_Median_AllSTN’ contains data for the median WithinR from each STN for the five models that are tested. (Data for Models 2 and 4 are presented in Figure 7).

https://doi.org/10.7554/eLife.19089.017

Download elife-19089-fig7-data1-v2.mat

When gripping was performed over a range of efforts we found that most information about the gripping force profile was contained in the beta and gamma bands in the STN LFP. The rich information in the gamma band was consistent with prior studies of grasp decoding using ECoG from human motor cortex (Flint et al., 2012, 2014; Pistohl et al., 2012). Changes in beta activity also contributed to the prediction of the force profile and inclusion of beta activity improved the encoding accuracy. This is consistent with our previous finding that beta desynchronization can encode gripping force, especially at the low effort levels (Tan et al., 2013, 2015). It is unlikely that the predictive power of activities in the beta and gamma band was related to contamination by movement artefacts as we used bipolar LFPs for the decoding in which any common artefact is removed through common mode rejection. In addition, we observed increased activity in the low frequency (delta and alpha bands), but features extracted from these frequency bands deteriorated the decoding of force during our core paradigm, even though separate parameters in a first order linear model were estimated for the low frequency activities. This contrasts with previous research which showed that the local motor potential and delta activity (0–4 Hz) in the LFP from the hand area of the primary motor cortex contains information about muscle activity (Flint et al., 2012) and pinching force (Flint et al., 2014).

Our first order linear dynamic model was cross validated in an independent patient cohort in a paradigm in which patients were asked to grip as fast and hard as they could. In this case, though, the additional incorporation of theta/alpha activity further improved the accuracy and reliability of the force prediction. This may relate to the focus on maximal voluntary contractions and their execution as fast as possible, placing additional attentional demands on the participants. Oscillatory activity in the theta/alpha range may be involved in attentional mechanisms (for review see Palva and Palva, 2007). In particular, alpha activity (7–13 Hz) in the subthalamic nucleus of patients with Parkinson’s disease is coherent with parieto-temporal cortical activity in a circuit that has been proposed to subserve attentional functions (Hirschmann et al., 2011; Litvak et al., 2011). Another consideration is the potential for the predictive power of low frequency activities to be related to contamination by movement artefacts secondary to movement overflow during maximal contractions.

Taken as a whole our results suggest that reciprocal changes in synchronised oscillatory activity in the STN can potentially provide control signals for the motor plant. In line with this is the impaired effort scaling in subjects lacking beta and gamma related motor reactivity in the STN LFP. We speculate that the motor plant may behave as a first order linear dynamical system translating a basal ganglia effort signal to force for a given effector. The present analyses allow quantitative assessment of the importance of such a simple transfer function model of basal ganglia-motor function in terms of its remarkable ability to predict both the pattern with which force is developed and the static force achieved on a single trial basis. Note that although the prediction is being used here in terms of its statistical meaning of accounting for the variance in a second signal, it also satisfies the physiological implications of this term in that changes in the STN LFP preceded changes in the measured force by about 200–300 ms. Our quantitative approach also allowed us to demonstrate that effort encoding is relatively lateralised in the basal ganglia.

Patients with Parkinson’s disease have been shown to have attenuated power modulation in both the beta band (Doyle et al., 2005; Devos and Defebvre, 2006; Androulidakis et al., 2007; Anzak et al., 2012) and the gamma band (Androulidakis et al., 2007) during movement initialization and when a constant force is meant to be sustained (Tan et al., 2013) when off dopaminergic medication. Such a diminished range of power modulation in the beta and gamma bands during movement may restrict the dynamic range of the coding of force in gripping in Parkinson’s disease. This may impair the force generation, and thereby cause bradykinesia, hypokinesia and weak movements when off compared to when on medication, presuming that the relationship between the STN force-encoding signal and generated force remains fixed across drug states (Figure 8). This is consistent with previous observations that untreated PD patients produce normal muscle activation patterns, but that muscle activity is not adequately scaled to produce the required force (Berardelli et al., 1986; Turner and Desmurget, 2010), PD patients are weaker off dopaminergic medication or DBS (Corcos et al., 1996; Alberts et al., 2004) and that PD patients show an increased probability of selecting slow movement speeds (Mazzoni et al., 2007). On the other hand, if the range of the maximal force to be achieved is to be kept similar to normal when off medication, the scale between the STN force-encoding signal and gripping force will be steeper when untreated (as indicated by Figure 8B). This will lead to abnormally high grip forces in Parkinsonian patients, as observed when they are asked to lift and hold an object (Fellows et al., 1998). Increased force scaling can also lead to the more coarse control of force and difficulty in finely tuning generated force. This is consistent with abnormally increased gripping force observed in precision grasping movements in Parkinson’s disease (Wenzelburger et al., 2002).

Figure 8

Download asset Open asset

Our findings suggest that the STN LFP could provide a high-performance control signal for BMI driven neuroprosthetic grasping in paralysed patients, leveraging advances in surgery for deep brain stimulation which has now become a relatively safe procedure (Larson, 2014). Surgical subcortical targets, such as the STN and globus pallidus (GPi), are involved in motor planning and execution. Activities from STN and GPi have been shown to correlate with movement parameters such as movement amplitude and speed (Brücke et al., 2012; Joundi et al., 2012), and are also modulated by movement intention (Kühn et al., 2006). Basal ganglia output has been theorized to regulate movement gain in healthy motor control, and can contain important information about the motor vigour (Shadmehr and Krakauer, 2008; Turner and Desmurget, 2010). Recordings of STN LFP signals have been shown to be stable over months using an implanted amplifier (Quinn et al., 2015; Neumann et al., 2016), and the signals are similar at re-operation several years later (Giannicola et al., 2012). It is also possible that the predictive potential of the STN LFP might be further improved in patients without Parkinson’s disease, as dopamine increases the reactivity of beta and gamma activities around the time of movement (Doyle et al., 2005; Androulidakis et al., 2007). Moreover the reliability of predictions was dependent on the relative scale of movement-related power changes in the gamma band which could be improved by refinements in electrode design and targeting, perhaps using this reactivity to select the optimum implantation target. Therefore, LFP signals recorded from the basal ganglia have the potential to provide stable and high-performance BMI input signals for the control of neuroprosthetic devices in paralysed patients, complementary or alternative to spikes or ECoG signals recorded from the cortex. In this regard, it is important to note that some LFP reactivity is retained when movements are imagined rather than actioned (Kühn et al., 2006), suggesting that peripheral afferents may not be necessary for the beta suppression and gamma increases tracked here.

Nevertheless, several critical issues need to be addressed before neuroprosthetic control by STN LFPs can be considered for implementation in BMIs. Any application is predicated on the assumption that the spectral reactivity demonstrated in the STN of treated PD patients is more-or-less preserved in chronically paralysed patients. This remains to be proven. Although our results suggest that basal ganglia LFP changes might potentially be useful in force control, the selection of movements or effectors is also important in BMI control. Whether and to which extent movements involving different body parts, such the lower as opposed to the upper limb, can be decoded from the STN LFP remains unaddressed, although microelectrode recordings suggest some spatial segregation in activity related to different limbs in the human STN (Rodriguez-Oroz et al., 2001; Theodosopoulos et al., 2003). In addition, we have not demonstrated the specificity of STN LFPs for the decoding of force. Finally, plasticity of cortical control signals is likely to play a significant role in the maximisation of the performance of BMIs (Carmena et al., 2003; Ganguly and Carmena, 2009), and yet whether subcortical signals can adapt over time remains to be seen.

There are some limitations in the current study which need to be acknowledged. First and foremost, we were unable to predict the force in about half of our patients. There may be many reasons for this, including disease related impairments despite dopaminergic therapy, post-operative stun effects and failure to pick up LFP activity from the ‘motor’ STN due to electrode targeting error. Confirmation that some electrode contacts were in or touching the STN was given by the surgical team at each centre upon review of pre- and post-operative imaging blinded to the electrophysiological results. According to this standard seven out of the nine electrodes affording force prediction were on target, with the other two electrodes in the posterolateral tail of STN. This suggests that significant movement-related spectral modulation in the beta and gamma bands, -the basis for selecting out this group in the first place, might be a good electrophysiological marker of proximity to the STN, and perhaps to the region of the STN involved in the motor representation of the hand in particular. Noteworthy in this regard, the spatial gradient of gamma reactivity was more acute than that of beta reactivity. In contrast, six out of the 12 electrodes that did not allow reasonable force prediction were identified as neither being in nor touching the STN and these subsequently revised on four sides or left inactive on two sides. A stun effect or disease related impairment might help account for the lack of joint beta and gamma reactivity in the remaining patients in this group, although two of these still had significant modulation in the beta band and one in the gamma band.

Such factors are also likely to have contributed the variation in prediction performance between subjects in whom some force prediction was possible. Amongst these, predicted force based on the STN LFP features was also noisier than the measured force. This could be caused by dynamic fluctuations and short-time-scale events in the STN LFP oscillatory activity (Feingold et al., 2015). It could be improved by incorporating a time history of the LFP signals which is equivalent to filtering the control input from the brain, and by incorporating filtering algorithms in the predicted force (the plant output of the model) in the time domain, such as the Kalman filter, or by convolving the output with a static nonlinearity (Fagg et al., 2009; Flint et al., 2014). Alternatively the relatively noisy predictions might in part reflect involuntary dyskinetic movement that did not impinge on grip force and yet may have been parameterised in the LFP control signal. Second, in a minority of trials grip forces were predicted based on the STN LFP but there was no measured force. This might represent a movement intention with consequent LFP change without movement execution, or voluntary movements that were not captured by the force dynamometer. Likewise these factors might have modified baseline LFPs in some trials leading to the negative force predicted in rare trials, unless predictions were clamped so as not to fall below zero. Third, the force prediction presented here focuses on the rest, force initialization, development and stable force holding phase of a grip, but does not include the force release or termination phase.

Finally, despite the observation that the first order linear model provides a good approximation of the relationship between STN LFP features and gripping force, the model may be an over simplification. Moreover, the results do not settle the discussion as to which circuits predominantly account for the selection of motor scaling, which can equally be attributed to cortical function (DeLong and Wichmann, 2010).

Despite these limitations, our findings do suggest that signals elaborated in and/or transmitted through the basal ganglia, and the STN in particular, carry information about motor scaling. We have also shown that features in the STN LFP combined with a simple dynamic model can be used to reliably predict the gripping force profile of the contralateral hand, even in individual grips. We propose that the LFPs from deep brain structures such as the STN could potentially provide stable and high-performance BMI input signals, complementary or alternative to neuronal spikes or ECoG signals recorded from the cortex. The recording stability and rich information content in the STN LFP about movement intentions and parameters make it an interesting signal with respect to BMI control.

In the main paradigm, eleven patients with idiopathic Parkinson’s Disease (mean disease duration 11.3 years, mean age 61.3 years, range 49–73 years; seven males) provided informed consent to take part in this study, which was approved by the local ethics committees. Patients underwent bilateral implantation of DBS electrodes into the STN, as a prelude to therapeutic high frequency stimulation for advanced idiopathic PD with motor fluctuations and/or dyskinesia. Techniques to target and implant electrodes in the STN have previously been described (Foltynie and Hariz, 2010). Microelectrode recordings were not made during surgery. The permanent quadripolar macroelectrode used was model 3389 (Medtronic Neurologic Division, Minneapolis, MN, USA) featuring four platinum-iridium cylindrical surfaces. Its contacts are numbered 0, 1, 2, and 3, with 0 being the most caudal and contact three being the most cranial. Localisation was supported intra-operatively by the effects of direct stimulation and by immediate post-operative stereotactic imaging. Nonetheless, in acknowledgement of the fact that not all electrode contacts could be expected to lie in the STN per se, we term the area sampled by the electrode contact the STN region (STNr). DBS electrode extension cables were externalized through the scalp to enable recordings prior to connection to a subcutaneous DBS pacemaker, implanted in a second operative procedure up to seven days later. One out of the eleven patients (case 2) had only one electrode externalised for testing, thus we could record from 21 STN regions (STNr). Clinical details of the patients are given in Table 1. The patients showed 53.4 ± 6.0% (p<0.001) improvement in the motor section of the Unified Parkinson’s Disease Rating Scale (UPDRS) on treatment with levodopa, indicating good responsiveness to this drug.

Subjects were seated in a comfortable chair with their shoulders adducted and their elbows flexed at about 90°. Subjects were first asked to grip the dynamometer with maximal effort three times, with each trial lasting for 3 s. Then they were presented with a series of imperative visual cues (red light-emitting diode illuminated for 3 s), separated by 11–13 s, and instructed to ‘choose an effort level from the scale provided and then to squeeze the force dynamometer at this chosen effort level when the light comes on and maintain this squeeze for the duration of the light’. Subjects were provided with the Rated Perceived Exertion Scale with 11 levels ranging from zero to 10 (Borg, 1998) printed on a piece of A4 paper. They were asked to try and randomise their selection of effort levels, so that the levels were varied from trial to trial, and all levels were represented. The subjects reported the effort level verbally after each grip. The mean number (± SEM) of trials per hand per subject was 31 ± 2 grips, with a mean number of trials per level per subject of 3 (±1). In particular, 2–3 trials were self-rated as maximal effort and 2–3 trials self-rated as minimal effort in each subject. Patients were asked to grip following illumination of the LED, but were not requested to respond as quickly as possible. There was no other feedback provided to the patients related to the generated force.

Recordings were made when the patients were ON their usual dopaminergic medication, 3–6 days postoperatively, while electrodes were externalized and before implantation of the pulse generator. Grip force was measured one hand at a time using an isometric dynamometer with standard Jamar design, and its handle set in the second of the five discrete grip diameter adjustments possible (G200; Biometrics Ltd, Gwent, UK). The order in which left and right hands were tested was counterbalanced across subjects. Monopolar LFPs recorded with a TMSi porti (TMS international, Netherlands) and its respective software. A common average reference was used for the monopolar recordings and these were low and high pass filtered at 0.5 and 500 Hz, respectively. Bipolar signals were derived offline by subtracting the monopolar recordings between neighbouring contacts on each electrode. The force was only a low pass filtered at 200 Hz. EMG signals from the first dorsal interosseous (FDI) of the activated hand were recorded in all patients, and the EMG signals from the extensor digitorum communis or extensor carpi ulnaris of both lower arms were recorded in six out of the 11 patients. LFP, EMG and force measurements were initially sampled at 2048 Hz. The effort level the subject reported verbally after each grip was logged manually and then used to label each individual trial.

Each electrode has four contact points, and the LFP data were converted off-line to give three bipolar contact pairs (01, 12 and 23) per electrode. Nonetheless, in acknowledgement of the fact that not all electrode contacts could be expected to lie in the STN per se, we term the area sampled by the electrode contact the STN region (STNr). Continuous wavelet transform, with Morlet wavelet and cycle number of 7, was then applied to LFP recordings from each bipolar contact pairs for time-frequency decomposition. The average power changes relative to the pre-movement baseline over the three trials of maximal effort gripping was calculated for each bipolar contact in the contralateral STNr. Three features from each bipolar LFP signal were extracted: the power change in the theta/alpha (4–12 Hz) band, beta (13–30 Hz) and gamma (55–90 Hz) frequency bands. For each electrode, the bipolar signal with the largest movement related reduction in the beta power was selected for analysis. Further, the average power in each frequency band and at each time point was compared against the distribution of the average power for that frequency over a one second period of time before the cue. Significant movement-related modulations were defined as those trials in which there were at least 50% of time points during the second after movement onset with power smaller than the 5% boundary (to capture event related desynchronization) or larger than the 95% boundary (to capture event related synchronisation) of the power distribution before movement for that frequency band. This procedure identified 9 STNs (from six different patients) with significant movement-related modulations in both the beta band and the gamma band; three more STNs showed significant movement-related modulations in either the beta or gamma band. However, there was no significant movement-related modulation in either beta or gamma band in the remaining 9 STNs. Details of the patient and the average power changes during the second after cue onset for maximal effort gripping in the beta and gamma bands for each STN are presented in Table 1.

For force decoding, time-frequency decomposition using continuous wavelet transform was applied to the STN LFPs from the bipolar contact previously selected. For each individual trial of gripping, the power change for each frequency at each time point was calculated by normalizing the power at that time point against the average power during the 1 s before cue presentation, so 0 indicates the power being the same as the baseline activity before cue, positive values indicate power increase (referred to as ERS) and negative value indicate power decrease (referred to as ERD). Then average power changes in theta/alpha ($\alpha$: 4–12 Hz) band, the beta ($\beta$: 13–30 Hz) and gamma ($\gamma$: 55–90 Hz) frequency bands at each time point were calculated. The latter frequency range of 55–90 Hz was selected on the basis of our previous study showing that STN LFP activities within this range increase during the onset of a grip (Anzak et al., 2012) and correlate with the stable force achieved during a grip (Tan et al., 2013).

Our previous study (Tan et al., 2013) showed that activity changes in the beta and gamma bands make major contributions to the encoding of efforts in gripping. In the present study, we tested different hypotheses about the dynamic relationship between LFP features and force. Based on the results from the previous study showing that the difference between gamma ($\gamma$) and beta ($\beta$) modulations correlates with effort at the force holding phase, the first model to be tested used the signal of $\gamma -\beta$ as the as a control input:

Where $\beta$ and $\gamma$ are beta and gamma band activity change in the STNr LFP, respectively. $\frac{K_p}{T_p\bullet s+1}{e}^{-T_d\bullet s}$ is a standard representation of a first order linear dynamic system with a time delay in the Laplace domain, where $K_p$ is the steady state proportional gain, $T_p$ is the time constant which is a measure of how fast the force output responds to brain signal changes, and $T_d$ is the time delay between the brain control signal (LFP) and the measured force which describes the latency between the timing of brain signal changes and force output onset, s= j ⋅ ω which is a complex variable. The equivalent of Model 1 in the time domain is: $Force\left(t\right)+T_p\bullet \frac{\partial Force\left(t\right)}{\partial t}= K_p\bullet (\gamma \left(t-T_d\right)- \beta \left(t-T_d\right))$, which suggests that beta and gamma modulations in STN LFP encode the instantaneous amplitude and the differentiation of force over time.

The second model to be tested assumed that the activities from beta and gamma bands have a different proportional gain, but have the same dynamic relationship in terms of time constant ($T_p$) and time delay ($T_d$) in encoding force:

The third model tested assumed that the dynamic relationship between force and activities at different frequency bands are different, with different values for the proportional gain ($K_p$), time constant ($T_p$) and time delay ($T_d$). Thus the generated force is the sum of the distinct processes with different inputs and different transfer function parameters (Model 3):

These models were compared against models which include extra information from the STNr LFP in the form of the relative power change in the theta/alpha frequency band ($\alpha$) (Model 4 and Model 5):

Models with only activities from single frequency bands were also evaluated to see if the combination of activities from different frequency bands was necessary for the decoding for force:

The predicted force was clamped so it did not fall below zero. The parameters in different models ($K_p$, $T_p$, $T_d$) were estimated for each STN separately.

Cross-validation was used to evaluate the generalisability and decoding performance of the proposed models with the STNr LFP features as inputs for predicting gripping force. For each STN/contralateral hand, the parameters of different models were first identified using least-squares optimisation applied on a set of training data to fit the corresponding STNr LFP features and measured force of the training data. The decoding accuracy of the models was evaluated by applying the model with identified parameters on another set of testing data. When evaluating the performance of the model applied to across-trial averages, the average of STNr LFP features and forces over multiple trials of low effort levels was first used as the training data to identify the model parameters, and the model was then tested on average data from high effort trials; or vice versa, i.e. using the across-trial average of high effort trials as training data and the averages from low effort trials as testing data. Bayesian information criterion (BIC) was used for model selection. The BIC value for each model is calculated as: $BIC=n\bullet \ln \left({\sigma }_e^2\right)+k\bullet \text{ln}\left(n\right)$; where ${\sigma }_e^2$ is the error variance: ${\displaystyle {\sigma }_e^2=\frac{1}{n}\cdot {\sum }_{t=1}^n{\left(F\left(t\right)-\hat{F}\left(t\right)\right)}^2}$ with $F\left(t\right)$ and ${\displaystyle \hat{F}\left(t\right)}$ the actually measured and predicted force of each time point respectively; n is the total number of time points; k is the number of free parameters in the model.

When evaluating the performance of the model on individual trials, a five-fold cross validation was used. All data from each recording session were partitioned into five equal folds. During each iteration, one fold was retained for testing, and the other four folds were used as training data to identify model parameters. Five iterations would allow for each observation being used for validation exactly once. The five results were then combined to produce a single complete estimation for one session. Different evaluation parameters including within trial correlation (WithinTrialR), correlation between the measured and predicted stable force were quantified based on the testing data.

In addition, for cross-subject validation, the 1^st order dynamic model with STN LFP features as inputs was used to predict the gripping force of individual trials in another independent patient group. Ten patients were recorded in the study when the patients were asked to grip as fast and as hard as they can in each trial in response to an external cue. The details of patient information and experimental paradigms were reported in a previous publication (Anzak et al., 2012). Twenty trials were collected for each hand when the patients were on their normal dopaminergic medication. Within each STN in this independent patient group, a 4-fold cross validation was used. For each iteration, 15 trials were used to fit the model to estimate model parameters and the model was used to predict force on the remaining five trials. The same procedure was repeated four times so each trial was used for prediction exactly once. The correlation coefficients between the predicted force and measured force for each individual trial were quantified.

All analyses were performed in Matlab (version 2012b). Median and the range of values are reported if the sample number is smaller than 10, or the distribution is not normal. Otherwise, means ± standard error of means (SEM) are presented throughout the text. Correlation coefficients were Fisher-z transformed before any statistical test, but the raw values were presented in the text and the figures.

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

This work was funded by the EU grant FP7-ICT-610391, Medical Research Council, Rosetrees Trust, and the National Institute of Health Research, Oxford Biomedical Research Centre.

Human subjects: Informed consent and consent to publish was obtained from patients before they took part in the study, which was approved by Oxfordshire Research Ethics Committee.

© 2016, Tan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.