Abstract
Objective. Deep brain stimulation (DBS) provides dramatic tremor relief in patients with severe essential tremor (ET). Typically, the VIM nucleus is the most effective brain area to target for high-frequency electrical stimulation in these patients. Correlation analysis between electrical local field potential (LFP) recordings from the thalamic DBS leads and electrical muscle activity from the contralateral tremulous limb has become an attractive practical tool to interpret the LFPs and their association with the tremulous clinical manifestations. Although functional connectivity analysis between brain electrical recordings and electromyographic (EMG) signals from the tremor has been of interest to an increasing number of engineering researchers, there is no well-accepted tailored framework to consistently characterise the association between thalamic electrical recordings and the tremorogenic EMG activity. Methods. This paper proposes a novel framework to address this challenge, including an estimation of the interaction strength using wavelet cross-spectrum and phase lag index while demonstrating the statistical significance of the findings. Results. Consistent results were estimated for single and multiple trials of consecutive or partially overlapping epochs of data. The latter approach reveals a substantial increase on the range of statistically significant dynamic low-frequency interrelationships while decreasing the dynamic range of high-frequency interactions. Conclusion. Results from both simulation and real data demonstrate the feasibility and robustness of the proposed framework. Significance. This study offers the proof of principle required to implement this methodology to uncover VIM thalamic LFP-EMG interactions for (i) better understanding of the pathophysiology of tremor; (ii) objective selection of the DBS electrode contacts with the highest strength of association with the tremorogenic EMG, a particularly useful feature for the implementation of novel multicontact directional leads in clinical practice; and (iii) future research on DBS closed-loop devices.
1. Introduction
Essential tremor (ET) is defined as a neurological disorder that causes involuntary abnormal repetitive shaking. This shaking can appear in different parts of the body such as the hands, forearms, or head [1]. ET is a common movement disorder affecting around four out of 100 adults over 40 years of age. It is considered to be a centrally driven tremor, and the constituents of the network comprise anatomical areas of the physiological motor system [2]. Treatment of ET is mainly based on pharmacotherapy and surgery for medically refractory cases [3]. The main surgical approach currently in use consists of continuous deep brain stimulation (DBS) through implantation of a depth electrode in the area of the ventral intermediate (VIM) nucleus of the thalamus, a key area in the neuronal loop generating the tremor [3, 4]. The depth electrodes after DBS surgery emit continuous electrical signals coming from the neurostimulator or pacemaker [5]. To confirm electrode position within the VIM thalamic nucleus, in addition to neuronavigation techniques to target this specific anatomical area, some centres perform microelectrode recordings during surgery to detect tremor-related electric neuronal bursts and find the ideal position for electrode implantation [6]. Through the routinely used macroelectrodes in DBS surgery, neuronal network electrical activity can be recorded from the VIM thalamus in the form of local field potentials (LFPs) a few days after surgery.
2. Material
In this study, we use recordings from two patients (aged 64 and 53, both female) who underwent DBS surgery for medically refractory ET. Ethics approval for use of patients’ EEGs and LFPs for the development of new quantitative EEG (qEEG) methods was obtained both from the University of Sheffield and the NHS ethics committees (SMBRER207 and 11/YH/0414). The multichannel Natus Quantum Amplifier (Optima Medical Ltd.) at a sampling rate of 16,384 Hz was used for all EEG/LFP/EMG polygraphy recordings (analogue bandwidth 0.01–4000 Hz).
In our institution (Royall Hallamshire Hospital, Sheffield Teaching Hospitals, NHS Foundation Trust), every patient undergoing DBS surgery for tremor is offered, five days after implantation of the depth electrodes, a comprehensive electrophysiological analysis to include LFPs from the VIM thalamus, electroencephalography (EEG), and electromyographic (EMG) polygraphy recordings, measuring the electrical activity from the muscles, to look at the correlation between the tremor and the thalamic network oscillations. The cortical scalp EEG recordings are not used in this work. There are four electrode contacts on each macroelectrode placed in the thalamus in close proximity to each other (contacts 0, 1, 2, and 3) through which bipolar recordings of LFPs can be obtained. Many centres use an empirical approach in selecting the ideal electrode contact to stimulate after DBS surgery to obtain the ideal tremor suppression. This approach can be time-consuming, and it will become more of a problem as new multicontact directional leads make their way into clinical practice [7]. Additionally, it seems that the LFPs offer very significant information which in the future could be used to optimise clinical outcomes in closed-loop systems [8]. However, before this becomes possible, an objective measure of the strength of association between the VIM thalamic network and the tremor recorded on EMG, or the equivalent mechanical oscillations through accelerometers, is required. One practical solution is based on the correlation analysis between data coming from the DBS leads (LFP) and electrical muscular activity (EMG) of the tremulous limb. Corticomuscular functional connectivity is defined as the interaction between the electrical activity of the cortex in the brain or the thalamus in this instance and the electrical activity recorded from various muscles. During the last few years, a few methods have been developed to understand this type of functional connectivity [9–11].
From the engineering perspective, the approaches to study connectivity between two signals can be classified as time- and frequency-domain-based methods. One of the classic linear measures to estimate similarity in time domain is cross-correlation that has been used to study ET [12]. It measures the degree of connectivity when a time series is shifted from other reference series. This method is also useful to measure the time lag between two signals. However, this method usually assumes that the system is linear and stationary. Another branch to quantify the correlation or causality between signals in time domain is by mutually predicting selected observable measurements based on multivariate modelling, where the best-established methods are based on the Granger causality test [13] that has been used to study the correlation between EEG signals [14]. They are based on parametric modelling. A full and unbiased model is therefore required, which can be a challenge to achieve due to the limited knowledge on the human brain. Other nonparametric methods include mutual information [15] and transfer entropy [16, 17], which are model-free but usually require larger datasets or averaging over many realisations to mitigate the effects of noise. In the scope of frequency-based methods, cross-spectrum allows determining the connection between two stationary signals in terms of frequency [18]. When computing complex cross-spectrum, results can be divided in cospectrum (in-phase connectivity) and quad-spectrum (out-of-phase connectivity). Coherence [19, 20] uses a normalisation of cross-spectrum values that takes unit value when total linear phase relationship is detected or zero value for independent signals. Both methods are easy-to-use and computationally inexpensive, but the investigated correlation must be stationary and linear and there is no-coupling information among various frequencies. Cross-bispectrum [21] is used to detect the quadratic phase coupling (QPC) between frequency components of two target signals. One pitfall of this algorithm is that bispectrum is affected not only by nonlinearity but also by non-Gaussian data. Thus, it requires the data to present Gaussian distribution to make sure to detect nonlinearity. A suitable approach to overcome these non-phase-coupling peaks is through cross-bicoherence [22]. However, this method only favours the strongly phase-coupled signals, the ones that show QPC interactions.
In addition to the frequency- and time-based approaches, efforts have been made with the wavelet domain to study connectivity of brain networks. This approach aims to address limitations of the above methods on tackling dynamic systems by providing time-resolving value with accurate locality. Meanwhile, it is a model-free (nonparametric) measure, which reduces the requirement of a priori knowledge of the underlying model. Wavelet coherence has attracted increased interest on studying brain-related disorders. Jeong et al. [23] proposed to use wavelet energy and wavelet coherence as EEG biomarkers to distinguish Parkinson’s disease and Alzheimer’s disease. A wavelet coherence-based clustering of EEG signals has been developed to estimate the brain connectivity in absence epileptic patients [24]. It has also been applied to studies on autism [25], traumatic brain injury [26], schizophrenia [27], and poor sleep quality [28]. However, there is very limited research focusing on its application on ET. Furthermore, for the approaches based on wavelet coherence to understand brain connectivity, there are no well-accepted complete frameworks to understand cerebromuscular connectivity systematically.
Addressing these challenges, this paper develops a new wavelet-based correlation analysis framework combined by the estimation of connectivity strength, significance test, and phase-delay characterisation. It aims to better understand cerebromuscular interactions in a structured manner. It should be noted that this framework has the prospect to be applied on other applications of connectivity analysis, such as EEG-EEG and EEG-EMG.
3. Methods
Model-free (nonparametric) measures are chosen in this paper to study the correlation between LFPs and EMG because this biological system is so complex that it would require a substantial number of parameters and computation time to build a satisfactory parametric model. Additionally, there is no well-accepted analytical model to start with. A common deficiency when applying biomedical signal-processing tools is the assumption of stationarity. A typical practice for the estimation of spectrum distribution is segmenting long records of data and averaging calculations over segments. The main disadvantage of this practice is the inability of having time-resolved values. There is therefore no time resolution, and the dynamic behaviour of neuronal interactions cannot be revealed. One solution to overcome this issue is the use of wavelets.
Wavelet transformation makes a decomposition of a time series into a frequency-time domain. It uses convolution of a mother wavelet and its scaled and shifted versions. Among all the possible mother wavelets, the Morlet wavelet provides a good balance between frequency and time resolution and has been widely used in EEG and EMG research [29]. This study is restricted to this wavelet. Equation (1) corresponds with a normalised Morlet wavelet where the dimensionless frequency is set as 6 and dimensionless time is denoted by . Considering an observed time-serial , in continuous wavelet transform (CWT), the mother wavelet is modified by varying the scale , as shown in (2), so that with a discretised time domain of time step .
3.1. Wavelet Coherence
Coherence is one of the most widely used methods for measuring linear interactions. It is based on the Pearson correlation coefficient used in statistics but in frequency and time domain. It measures the mean resultant vector length (or consistency) of the cross-spectral density between two signals. Its squared value varies from 0 to 1, meaning low and high linear frequential correlation. During this study, coherence is used as a reference standard for comparison to other methods. The wavelet formulation of coherence between two signals, and , and in the frequency and time domain, can be formulated as where is the wavelet cross-spectrum between and and and are the corresponding autospectrums. Working with two single signals (single realisation) usually requires using a smoothing operator (see operator in (4)), and ergodicity properties should be assumed [30].
If multiple trials of both signals are available, square coherence can be estimated using (5), which is used in this study. where
The number of trails is denoted by .
To assert significant values of wavelet coherence, the statistical methodology stablished by Gallego et al. [31] is employed. It estimates a threshold based on the null hypothesis () of independency by analytically calculating the statistical distribution of coherence. Specifically, assumes that both signals are independent Gaussian variables. Under , with a specified probability , where is calculated as
In this paper, the parameter is set as a fixed value of , equivalent to a 95% of confidence interval.
3.2. Wavelet Cross-Spectrum
Wavelet cross-spectrum (WCS) also provides information about linear synchronisation, but its values are not normalised as in wavelet coherence. Its calculation is written in (6) for multiple trials. Its values seem difficult to be interpreted by statisticians and complementary plots, such as coherence or coquadrature, that can help understand the frequential relationships between signals [9]. That is why few neurological connectivity studies have used this method. However, this paper proposes to use an appropriate significant test that allows WCS results to be judged with more clarity by including the variance of each signal.
It has been proven by Bigot et al. [32] that including the variance spectrum data (autospectrum) and benefits the interpretation of WCS in contrast to wavelet coherence. Its performance is based, analogously to coherence, on a threshold calculation under a null hypothesis. However, this makes the statistical procedure using a combination of parametric and nonparametric estimations. It has parametric characteristics since it stablishes the signals as independent Gaussian vectors. But its distribution is considered to show a zero-mean value and a general covariance matrix and estimated by the sampled data. Hence, it does not make any parametric assumption on the covariance. The statistical test stablishes, under and , that . The threshold , for each time and frequency point, is defined as where and are the largest eigenvalues of the empirical covariance matrix of the time series x and y, respectively. The symbol means the energy of the wavelet, but since wavelet-normalised values are considered, this variable is omitted. This paper considers a significant probability of 95% with .
3.3. Phase Lag Characterisation
Apart from the quantification of the linear interaction strength and associated significance, it is also important to measure the time or phase lag between signals at a certain frequency. It is particularly important for studying corticomuscular interactions as they carry significant time delays.
The most straightforward approach would be calculating the time lag from the phase information in the complex data of the WCS results with a simple fraction where is the phase difference between both signals at a specific time and frequency and is the sample rate. However, volume conduction can cause the coherence and the phase-locking value to spuriously increase, and (9) is not effective to deal with this kind of common noise sources. To overcome this problem, a measure called the weighted phase lag index (WPLI) algorithm [33] that describes the consistency of the phase difference or time lag is used. This value will inform if the phase lag existent between both signals is consistent. This method not only highlights frequency bands where the phase difference is constant, indicating a strong linear relationship, but also penalises those synchronisations where time lag is close to zero, thus avoiding possible artefacts coming from common noise sources. The weight can be illustrated by Figure 1, which clearly shows the WPLI weights cross-spectra according to the magnitude of the imaginary component of the cross-spectrum. Cross-spectra around the real axis contribute to a less extent than cross-spectra around the imaginary axis. Their values are calculated as where is the imaginary part of the argument and is the cross-spectrum value of the th trial of the total number of trials (). The value is normalised from −1 to 1. As indicated by Figure 1, a value of 1 suggests that the phase lag is 90° while a value of 0 corresponds to 0° or 180°. It should be noted that this process should be applied before the significance test.
The proposed methodology can be illustrated by Figure 2. Starting from the multiple trials’ data collection or extraction, the WCS is applied on each trial to measure the correlation strength. After the number of trails is sufficient, the results of WCS of all trails are averaged. A significance test can then be applied to produce the binary correlation map. Meanwhile, the phase lag can be estimated based on the averaged WCS.
4. Results on Simulation Examples
The simulation example aims to evaluate the performance of the proposed analysis framework on frequential and time resolution and robustness against noise. Considering a linear single input single output (SISO) system, the input signal is defined as
The input signal constitutes an ensemble of components, with different frequencies in the cosine form. It is defined within a temporal range through a window function . The entire ensemble was sampled by a normal distribution . In other words, the amplitude of the entire signal is modulated by a random Gaussian signal. For the window function, a Tukey window was used with a tapper parameter of 0.5. This type of window was selected to avoid possible ringing artefacts coming from the temporal sharp transitions of the frequency components. A white Gaussian noise was also added, emulating a more realistic situation. The output signal is defined as
One difference with the input is the change on the amplitude of each frequential component by setting . In addition, each component presents a specific phase . Therefore, the linear system modulates the amplitude for each frequency of the input and delays it with a specific value. It must be noted that one can get several realisations of both measures, having then multiple trials of this linear interaction. Figure 3 illustrates the produced input and output signal of a single trial, which includes two frequential components: within and within . The total data length is 1 sec with the sample rate of 1000 Hz. Both amplitude gains are equal and unit valued (), and delays and are set to 25 ms and 14 ms, respectively. Gaussian noise was added for both signals with the signal-to-noise ratio (SNR) of 40 dB. The bottom graph of Figure 3 shows the result of standard coherence based on Fourier transform. Although it successfully reveals the strong coherence at the frequency around 10 Hz and 30 Hz, it cannot localise when the coherence changes occurred. Such an approach therefore misses the time-resolved information, which is key to study a dynamical system. Although such a limitation can be partly addressed through using a sliding-window technique [34], the selection of window size is usually challenging and depends on the frequency of signals.
If there is no noise, the measured frequency-time interaction after the significance test using WCS can be represented by Figure 4. The yellow colours, indicating the significant interactions, clearly correctly capture the two frequential components and corresponding starting and ending times, which cannot be revealed by the standard coherence. The white dashed line marks the cone of influence (COI). When computing CWT using convolution procedure, edge artefacts cannot be ignored. That is why COI is introduced to describe the area in which the power of the shifted wavelet drops to of the value at the edge [35].
If there is noise involved, the number of trials is important and should be considered. Figure 5 shows the result of WC with the associated significance test where SNR = −5 dB and . The linear synchronisation around frequencies of 10 Hz and 30 Hz is distinguishable. The temporal resolution of the interaction is close to the ideal one since the ending points of high coherent frequencies at 10 Hz and 30 Hz are located close to [0.1 s, 0.5 s] and [0.6 s, 0.9 s]. However, high coherence values can be observed in the high-frequency band and very low-frequency band, which are determined as significant, which are artefacts. This is caused by the introduction of severe noise. Figure 6 shows the results using WCS for the same parameter settings. It shows a clearer plot with better contrast comparing to Figure 5. The simulated synchronisation can be easily differentiated from the rest of the noisy values based on the graph of significance test (see Figure 6(b)). This is explained by the normalisation process followed in the CWT process by Grinsted, where the energy of the wavelet signals is normalised ( normalisation). One direct consequence of normalisation is that the values at higher frequencies are compressed more than those at lower frequencies.
(a)
(b)
(a)
(b)
To quantitatively study the influence of SNR and the number of trials on the results, the 2-D correlation coefficient between the significance tests of the ideal case (see Figure 4) and the testing case (see Figure 6(b)) is employed to measure the accuracy of connectivity detection. The equation to calculate the correlation coefficient of two images A and B can be written as where and denote the means of A and B, respectively, and and denote the indexes of horizontal and vertical directions, respectively. A closer value of correlation coefficient to 1 indicates a better performance of detection. In this test, the SNR was varied from 30 dB to −10 dB with the step of 1 dB, and n was varied from 1 to 50 with the step of 1. The result is illustrated by a coloured matrix of correlation values, as shown in Figure 7. Drawing a horizontal line across the SNR axis, it can be observed that the correlation value increases with the number of trials. The increment rate is faster at higher levels of noise (lower value of SNR) than lower noise levels (higher value of SNR). A distinctive high noise range exists from 5 dB to −5 dB where limited trials (less than 6) can produce good results. For SNR < −5 dB and SNR > 10 dB, a significant number of trials are required to achieve reliable results. The fact of getting better results with high noise levels (from 5 dB to −5 dB) with limited trials than those with low noise levels (>10 dB) can be explained by the influence of the Tukey window. When applied to short-time segments of low-frequency components, such a window function generates low-frequency artefacts. These artefacts can pass the statistical test when the threshold is low, which happens if the variance of both signals is low, or, in other words, when low noise is applied to the signal. This means, the noise level can influence the variance and thus the statistical threshold and finds a balance where the threshold only shows the sequential linear interactions without artefacts.
To complement the WCS results, Figure 8 shows the estimated WPLI results with the noise level of −5 dB and 5 trials. Although the first interaction at 10 Hz is still differentiable, the second interaction at 30 Hz is not so clear when comparing with the rest of the noisy values. It should be noted that WPLI values are not aiming at locating strong linear interactions but indicating the consistency of time lag, which helps distinguish the true interaction and artefact. Combining Figure 8 with Figure 6, it can be observed that the synchrony at 10 Hz shows higher phase lag consistency than that at 30 Hz. Considering the time delay of the two components ( and ) were set to 25 ms and 14 ms, respectively, the phase difference is and . According to Figure 1, the weight of is higher than the weight of ; hence, WPLI values around the interaction of have higher contrast than those of .
5. Results on Essential Tremor Data
5.1. Data Collection
All electrophysiological recordings were obtained with a multichannel Natus Quantum Amplifier (Optima Medical Ltd.). Four types of data were available for each recording: scalp electroencephalography (EEG), intracranial thalamic (VIM) local field potentials (LFPs), electromyography (EMG) with surface electrodes, and monoaxial accelerometer recordings from the hands and head. All data were sampled at 16.38 kHz and then downsampled to 2 kHz. This study focuses only on the interactions between thalamic LFPs and contralateral EMG (as each right and left half of the brain supplies the contralateral side of the body). LFPs refer to the summated electrical neuronal activity recorded with the DBS leads from the VIM thalamus. This activity was recorded by a quadripolar DBS lead with three possible input channels (0-1, 0–2, and 0–3) taking the pole 0 as a reference. Figure 9 shows an illustration of the lead layout. EMG data was recorded by surface EMG electrodes placed in different arm muscles. Given three LFP channels and five EMG channels, there are 15 possible pairs to be analysed in terms of interactions. However, based on previous observations from our clinical work, the right triceps brachii muscle (commonly showing well-formed tremorogenic oscillations) and the left 0–3 LFPs were considered in this work. Therefore, this pair is primarily studied below.
5.2. Multiple-Trial Data Extraction
The simulation example has demonstrated that the use of multiple trials or realisations of an event was important to better reveal significant interactions. However, it is not possible to repeat the same ET event in the same subject with the same conditions such as timing within each tremorogenic oscillation. To overcome this issue, the tremorogenic EMG signal trials were substituted by time segments coming from a long single trial signal. The necessary condition for extracting statistical properties by analysing data over time instead of evaluating several data samples is called ergodicity [36]. This is an important assumption that the target signal is considered as stationary instead of nonstationary and thus not assuming dynamical changes along the entire signal. However, it does not impede the detection of dynamic changes on frequential interactions. If the dynamic changes occur periodically, the true interactions will be enhanced by considering those time segments, while noise will be attenuated. Additionally, restricting time segments around a specific time interval, with the right number of trials and overlap, will make possible to see how local stationarity properties change over time, similar to a sliding-window technique, but the selection of widow size is not required.
For correcting this issue of the timing offset, a reference point, for each time segment extracted, is set as the closest peak of the EMG signal. Through this manner, phase cancellation artefacts can be avoided to some extent. To detect the peaks on EMG, the procedure follows three steps. The first step processes the data with a linear low-pass filter (passband edge frequency 15 Hz, stopband frequency 30 Hz, passband ripple 1 dB, and 60 dB of attenuation) since it is known that the tremor appears at low frequencies and the filtered signal is corrected with the corresponding delay of the filter. In the second step, the frequency component with highest magnitude is analysed. In the third step, a peak neighbourhood search is performed with a restriction based on the period of the fundamental frequency. Figure 10 illustrates an example of peak extraction from an epoch of EMG sample data coming from the triceps brachii, where the top figure shows a time series of raw EMG data and the bottom figure shows the filtered data with the detected peaks marked by green lines.
5.3. Interaction Estimation Based on a Single Trial
Figure 11 shows the result of WCS of a 10 sec single trial between left L0L3(LFP) and right triceps (EMG). As shown in the top graph, there is a distinguishable linear interaction between LFP and EMG around 5 Hz, which at the same time is the frequency of the arm tremor (validated by the accelerometer data). Moreover, the intensity of 5 Hz interactions is intermittent, showing peaks and troughs over time. In addition, weaker synchronisation is present at lower and higher frequencies. The significance test of a single trial is shown in the bottom graph of Figure 11, where significant interactions are observed from 1 to 128 Hz when using a single trial that as a result reduces the confidence level of estimation. Multiple trials are required to improve the confidence level to reveal the true significant interactions by reducing the influence of noise.
(a)
(b)
The colour scale establishes a range of values, in terms of WCS module, of fourth order of magnitude. However, it is not a reliable feature to compare with other combinations of signals, since the scale depends on the amplitude of the original signals that can differ depending on the impedance and position of the reference electrode. Instead, comparing different plots during a long period of time and checking which one is more consistent and regular will lead to better interpretation of the underlying interactions. Figure 12 shows the WCS results of a 15 sec single trial of five combinations between LFPs and EMG recordings. In this figure, the triceps brachii (a) presents a more stable and periodical behaviour in terms of the linear spectral correlation with L0L3. Figures 12(b), 12(c), and 12(e) show high spread values over lower frequencies.
(a)
(b)
(c)
(d)
(e)
5.4. Interaction Estimation Based on Multitrial
Two parameters to be determined for multitrial analysis for real data include the number of trials and the overlap rate. With sampling starting from the 6th sec of the data of L0L3 and right triceps brachii with a window length of 1 sec epoch, the correlation was estimated where the trial number was set as 10 and 20, and the overlap rate was set as 0%, 50%, and 75%. For the trial number of 10 and the overlap rate of 0%, the sampling windows are [6 s, 7 s], [7 s, 8 s], … , [15 s, 16 s]. It should be noted that the overlap rates of 50% and 75% are approximate values. The true overlap rates are determined by the references points based on the closest peak of the EMG signal. For example, the second window of the overlap rate of 50% is not necessary to start exactly from 6.5 s. It starts from the closest EMG peak around 6.5 s. Figures 13–14 show the results of WCS with the trial number of 10 and 20 and overlap rates of 0%, 50%, and 75%. In comparison with the result of a single trial shown in Figure 15, results of the multitrial highlight those interactions that keep repeating over time and reduce the occasional high values that could be artefactual. For example, the 5 Hz tremor along with its first harmonic at 10 Hz shows a stable linear correlation. Despite losing time resolution due to averaging, dynamic changes in WCS values at higher frequencies are still observable and could be genuine. A more regular temporal pattern can be observed for the frequency band from 16 Hz to 64 Hz following the increment of the overlap rate. The increased number of trials makes the results smoother and more consistent.
(a)
(b)
To evaluate the time lag of the significant interactions, the WPLI approach was applied and the result is shown in Figure 16. The higher contrast shows that components around the first harmonic of the tremor have the more prominent phase lag.
6. Conclusions
This paper proposes a novel data analysis framework to study thalamomuscular associations in essential tremor involving three steps: correlation strength estimation, significance test, and phase lag characterisation. This framework aims to improve the robustness and reliability of correlation analysis between the local field potential recordings from the brain and the tremulous electrical activity recorded on EMG. It has been shown in the simulation example that the proposed approach can effectively evaluate the linear interaction between two signals. The sensitivity analysis studies show how the number of trials and noise level of measurement affect the results. For data with noise level < −5 dB or >10 dB, a significant number of trials produce much better results. However, for data with noise level > −5 dB and <5 dB, the number of trials has less influence on the findings. The application of the method, on real data from two patients with ET undergoing DBS surgery for tremor suppression, demonstrates the validity of the proposed approach through segmenting a long single epoch into a number of overlapped windows to produce the averaged strength of associations. One limitation of this approach is that the result is difficult to be quantified due to the complexity of WCS patterns if the ground truth is unknown. Another potential limitation is that the number of trials plays an important role in improving the performance of this approach. With a single trial, the significance test cannot be constructed. In the real data application, it is not possible to repeat the same ET event (i.e., tremorogenic oscillation) in the same subject with the same conditions. Future work therefore will focus on quantification of the results and reduce the dependency from the number of trials.
It should be noted that wavelet cross-spectrum and phase lag characterisation used in this framework are not novel. However, combining them together along with a significance test is new. Furthermore, this is the first attempt to apply wavelet-based correlation analysis on patients with medically refractory essential tremor undergoing surgery. This paper shows a clear association between the thalamic local field potential recordings and the contralateral tremorogenic EMG oscillations, at the frequency of the tremor and its first harmonic (Figure 13). These interactions are beyond the observational empirical interpretation of thalamic LFPs and EMG recordings from patients with ET. This paper offers a framework that can be used to choose the thalamic contacts that show the strongest association with the tremulous EMG oscillations that if selected for high-frequency DBS stimulation produce a better clinical outcome in comparison to the empirical programming of the DBS device. Our clinical experience so far confirms this hypothesis, but this will be explored in a future study, on a significant number of patients with ET, where the emphasis will be on the use of such a framework to obtain best tremor suppression by objective selection of the ideal contacts to stimulate. As more complex multicontact directional DBS leads are already used in various centres, evidence-based (both radiological and neurophysiological) selection of the ideal contacts to stimulate will become increasingly important [7]. With this framework, the electrophysiological LFP recordings and their relationship to the tremor can be analysed to determine the DBS lead contacts that show the strongest association with the tremor and select them for stimulation. This work offers the proof of principle required to assess the utility and the limitations of this methodology. It has been demonstrated that the proposed framework can reveal significant cerebromuscular interactions, in this instance thalamic (VIM) LFPs versus the tremulous EMG activity, reaffirming in vivo that this part of the thalamus is part of the central tremorogenic network in ET. It could play an important role for future research on developing a closed-loop DBS device. It also has the potential to objectively determine in individual patients which of the thalamic lead contacts shows objectively the strongest association with the tremor, particularly as multicontact leads make their way into clinical practice. This thalamic lead contact could plausibly offer the best tremor suppression for each patient. This hypothesis will have to be confirmed in future electroclinical studies as our clinical experience so far is pointing in this direction.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the Through-Life Engineering Services Centre at Cranfield University (UK), in part by Royal Hallamshire Hospital at Sheffield (UK), in part by Zhejiang University of Technology (China), and in part by the Computational and Software Techniques in Engineering MSc Course at Cranfield University (UK). The authors thank Neurocare for the purchase of the medical equipment used for the polygraphy recordings on their patients. This is research was carried out in part at the National Institute for Health Research (NIHR) Sheffield Biomedical Research Centre (Translational Neuroscience)/NIHR Sheffield Clinical Research Facility. This work was also supported by the National Science Foundation Program of China (61601029), a grant of the Ningbo 3315 Innovation Team, and the China Association for Science and Technology (2016QNRC001).