The 1-3-3-1 leads, shown in Figure 1, are known as directional leads (see text footnote 2). Four levels (0, 1, 2, and 3) of electrodes are arranged along the length of a lead. All the available measurement directions are grouped into channels. There are 15 channels per lead, for a total of 30 channels if leads are implanted in both hemispheres (as is the case for the recordings in this paper). As explained in Figure 1, these leads allow stimulation and recording between levels (parallel to the leads) and between segments of a multi-segment level (lateral direction to the leads). The recordings were made with stimulation turned off. We download the LFP time series as a JSON file from the tablet that is used to communicate with the neurostimulator. The time series for each of the 30 channels consists of 21-s LFP recordings sampled at 250 Hz. Medtronic also applies two 100 Hz low-pass and two 1 Hz high-pass filters to the signal. We use these LFP signals for all the analyses in this paper.

Here we illustrate the phenomena of beating or amplitude modulation by using simple mathematical reasoning. Consider the superposition of two sinusoids with frequencies f₁ and f₂ (which lie close to each other) by using the identity,

which shows that the superposition of two waves with similar frequencies f₁ and f₂ produces a carrier wave that oscillates with the frequency,

and whose amplitude is modulated with a lower frequency wave with frequency,

The beat frequency is defined as the frequency of modulation of the power in the signal (square of the magnitude of the signal),

We note that the beat frequency f_beat is much lower than the frequencies f₁ and f₂, particularly when these frequencies are close to each other.

In this paper, Python v3.2.10 and its scientific packages are used for all signal processing and analysis of the signals. Fast Fourier Transforms (FFT) algorithm (McClellan et al., 2016) evaluated by using the scipy.fft.rfft function, is used to calculate the frequency spectra of the LFP signals as well as the limb tremor spectra obtained from the video of the patient. Frequency-bin sizes of 0.047 Hz for LFPs, and 0.05 Hz for the video recordings are employed. The square of the magnitude spectrum from the FFT only shows power at a particular frequency without any information about when a frequency event is happening in time. To get that information we also need the phase spectrum of the FFT. However, it is not possible to visualize time events by just looking at the magnitude and phase FFT spectra together. To visualize beta power bursts in detail we need to show the signal’s power both in the frequency and the time domains simultaneously while also obeying the time-frequency uncertainty principle (Addison, 2020; Mallat, 2008). We use the Continuous Wavelet Transform (CWT) scalogram of the LFPs time series for this visualization of the LFP signal and also the measured tremor time series. To concentrate on the beta band, the LFP signal is passed through a 13–31 Hz 5th-order discrete Butterworth band-pass filter. Butterworth filters are also known as maximally flat magnitude filters since they are designed to give a frequency response as flat as possible in the passband. Python’s function named cw with gaus8 wavelet from the pywt package was used to obtain the scalograms. We use 8 and 50 Hz for the upper and lower frequencies of the LFP scalograms; for the measured limb tremor scalograms, we use 1 and 10 Hz.

To make the detection of the dominant twin peak in the LFP spectrum objective, we use software to detect the peaks in the spectrum automatically with a desired prominence (quality) from all the 30 channels. We use the find_peaks function of Python’s scipy.signal package with the prominence parameter set to 0.04. Before detecting the peaks programmatically, the FFT of the signal needs to be smoothed appropriately to reveal the dominant peaks. Note that we are smoothing the FFT spectrum in the frequency domain to specifically observe and detect the peaks in the FFT. A low pass filter that is designed to give a flat passband when used on a time series, such as the Butterworth filter, is not appropriate for this purpose since it potentially destroys the spectrum’s tenacity and smears the peaks, making it difficult to use software tools to detect the peaks. Instead, we use the Savitzky-Golay filter popularized by Whittaker and Robinson (1924) and Guest (2012) for smoothing the spectrum. This filter, popular in the signal processing community, works by fitting subsets (windows) of adjacent data points with low-degree (order) polynomials by using the least square fit method³ (Riordon et al., 2000). A window size of 51 and polynomial order of 3 is used in this paper. When we say that we have detected a twin peak in the spectrum of a channel, it is implied that our software has detected twin peaks with a permanence equal to or greater than 0.04.

Google’s MediaPipe (see text footnote 1) is a framework based on AI for developers to build multimodal (e.g., video or audio) machine learning pipelines. Its primary use case is to provide machine learning models that can run live on smartphones or other edge devices. It can be used to rapidly analyze videos of patients taken by mobile phones or tablets. In this paper, we make use of MediaPipe’s “Pose Landmark Detection” (see text footnote 1). It uses its internal AI models to automatically identify key location points on a human body in a video or an image. These points are located at major connections to limbs and torso, such as wrists, heels, and knees. The possible location points are illustrated in Figure 2 on a human body (see text footnote 1). The framework also automatically evaluates a time series of the x and y coordinates of the key location points from a video. This information is sufficient to measure the patient’s tremor frequency, f_{tremor_measured}; either the x or y coordinate time-series can be used, and both give the same frequency characteristics. Since we are not interested in measuring the amplitude of the tremor in this paper, we do not require a time series for the z coordinates to account for the total motion trajectory of the limb. The limb tremor measurements are all obtained from a patient’s video made on December 15, 2022. A Sony HDR-CX440 Handycam video camera mounted on a tripod is used. The patient is seated on a sturdy chair with no other sources of vibration. The framework is used on this 854 by 480-pixel video, sampled at 30 fps (frames per second) to determine the x and y coordinates of the right and left wrists, knees, and heels of the PD patient.

In this paper, we use 21.15 s of channel 17’s LFP recorded signal with 5288 time-samples. The sampling rate is 250 Hz. Figure 3 shows the magnitude spectrum of this signal obtained by using the Fast Fourier Transforms (FFT) algorithm (McClellan et al., 2016). The bounds of the alpha, beta, and gamma spectral bands are shown by red dashed lines. For the boundaries of the spectral bands, we use 8 Hz, 13 Hz, and 31 Hz for the alpha, beta, and gamma bands, respectively. Note that it is hard to detect dominant peaks in the signal, especially by software tools, because of the lack of smoothness.

Having laid the theoretical background, we now address the central problem mentioned at the beginning of this paper and connect the LFP spectrum to the spectrum of tremor by formally forming a hypothesis in three parts.

(a) We hypothesize that the frequency of the power burst trains, in the LFP scalogram (which we have shown to be f₁ − f₂, is equal to the tremor frequency f_{tremor_measured}. In fact, the limb tremor of a patient is not composed of a single frequency, instead, the tremor spectrum consists of a peak located at f_{tremor_measured} with a non-zero width Δf_{tremor_measured}. Therefore, we calculate f₁ − f₂ over all the channels exhibiting dominant twin peaks to evaluate the mean (f₁−f₂)_mean. The hypothesis associates this mean value with f_{tremor_measured},

(b) The second part of the hypothesis links Δf_{tremor_measured} to the LFP spectrum. We evaluate f₁ − f₂ over all the channels exhibiting dominant twin peaks and calculate the standard deviation (f₁−f₂)_{std dev}. The hypothesis associates the standard deviation with Δf_{tremor_measured},

(c) We show later that the tremor fluctuation time, T_{tremor_fluctation}, can be evaluated from Δf_{tremor_measured}. The third part of the hypothesis states that T_{tremor_fluctation}is equal to the time duration, T_{lfp_groups}, of groups of beta-band trains seen in the scalograms of the LFP (top plot, Figure 4),

In summary, our hypothesis connects (predicts) three defining features of the measured tremor spectrum, f_{tremor_measured}, Δf_{tremor_measured}, and T_{tremor_fluctation}, to three independent features, (f₁−f₂)_mean, (f₁−f₂)_{std dev} and T_{lfp_groups}, seen in the LFP signals.

In this section, we investigate the distribution of energy of the LFP spectrum among the alpha, beta, and gamma bands for all the 30 channels. We define the LFP energy in a band as the integrated power of the spectrum over the frequency range of that band. Figure 11 (top plot) shows the distribution of the energy for all the 15 channels of the left lead, and Figure 11 (bottom plot) shows the same results for the 15 channels of the right lead. The histograms are drawn to the same scale for all the 30 channels. For our patient, the LFP energy in the bands is much higher in the left channels than in the right channels.

Now we evaluate the mean (f₁−f₂)_mean, and the standard deviation (f₁−f₂)_{std dev} by using all the eight channels that show Dominant Twin Peak activity (marked by TP in Figure 11). The results are shown in Figure 12, with a mean (f₁−f₂)_mean of 5.26 Hz, and a standard deviation of deviation (f₁−f₂)_{std dev} of 0.31 Hz. Hence, based on Equations 5a, b, our hypothesis predicts a tremor frequency for the PD patient to be 5.26 ± 0.31 Hz. We note that this width, 0.31 Hz, has no obvious relationship with the width of the peaks in the beta band of the LFP spectrum.

Google’s MediaPipe (see text footnote 1) framework was used to obtain a time series consisting of 550 time-samples (18.33 s) of the x and y coordinates of the limbs from a video of the patient. Figure 13 shows the tremor power spectrum (square of magnitude FFT) obtained from the time series of the PD patient’s right hand’s wrist movement. The spectrum shows a right wrist tremor peak located at f_{tremor_measured} = 5.22 Hz, and a Half Width at Half Max (HWHM) width Δf_{tremor_measured} = 0.32 Hz. Predictions made by our hypothesis, Equations 5a, b, are shown in green. Note that the measured quantities f_{tremor_measured} = 5.22 Hz, and Δf_{tremor_measured} = 0.32 Hz, show a close match with the predictions of our hypothesis 5.26 ± 0.31 Hz obtained from the LFP (Figure 12).

We restrict ourselves to the beta-band and ask this question: Can the spectral twin peaks observed in this paper be caused by the presence of random noise in the spectrum? In other words, if we start with a magnitude spectrum which contains one dominant peak at frequency f₁, with prominence ≥ 0.04, then what is the statistical likelihood that noise in the spectrum alone can cause another dominant peak with frequency f₂ to appear with prominence ≥ 0.04, which also satisfies the condition that |f₁−f₂| lies within the range:

Where σ_f (standard deviation of the measured peak) was evaluated by using the relationship HWHM = Δf_{tremor_measured} = 1.2σ_f mentioned earlier in the last section, and Half Width at Half Max (HWHM) was read off from the measured magnitude spectrum (Figure 13).

We address this question in two different ways by performing a pair of simulations comprising a large number of surrogate trials. In each of these trials, we simulate the noise in the spectrum to have a Root Mean Square (RMS) value consistent with noise exhibited in our LFP spectra; specifically, we choose channel 5 because it exhibits a dominant single peak with prominence ≥ 0.04. We calculate the LFP noise in the beta band of channel 5 by subtracting the raw noisy LFP spectrum obtained from channel 5 from the smoothed (Savitzky-Golay filter) LFP spectrum. The calculated RMS = 0.08 value of this subtracted signal (in the beta-band) is then used to generate uniform noise for all the trials in our simulations. As explained below, for simulatin-1 we superimpose noise onto an A/f^α spectrum (known as 1/f decay) fitted to channel 5, while for simulation-2, we add noise directly to the smoothed spectrum of channel 5 with the pre-existing dominant peak.

Simulation-1: We note that neural LFP signals commonly exhibit a 1/f decay in their power spectrum (Bédard and Destexhe, 2009), which follows a straight line on a log-log plot of the power spectrum vs. frequency. On a linear plot, this leads to a magnitude spectrum with a A/f^α (Smith, 2011, Chapter 1.7.3) behavior, where A and α are obtained by fitting to LFP data. We determine these parameters by fitting a straight line to the log(power spectrum) vs. log(frequency) plot of channel 5’s spectrum. This fit gives us a slope m = −0.85 and an intercept b = −1.21. From these, we obtain the parameters for the magnitude spectrum, A=e^b/2 = 0.55 and α = −m/2=0.42. Using a Python program, we generate 100,000 surrogate trials. A magnitude spectrum is constructed in each trial by generating uniform noise consistent with the desired RMS value and superimposing it onto the magnitude spectrum A/f^α. To detect peaks in each of the trials, we apply the Savitzky-Golay filter (as is done throughout this paper) and use our peak detection algorithm in the beta band with the prominence parameter set to 0.04. To detect twin peaks, we use the criterion |f₁−f₂|≤[4.68, 5.76]Hz. Out of a total of 100,000 surrogate trials, we extract the following statistics (CI stands for 99% Clopper–Pearson Confidence Intervals (CI) based on the 100,000-trial simulation):

512 trails are detected with twin peaks among 56,189 trails with at least one peak with prominence ≥ 0.04.

Conditional probability of observing a twin peak given that at least one prominent peak exists in the magnitude spectrum: 0.009112 (0.91%)

99% CI: [0.00811, 0.01019]

512 trials are detected with twin peaks among a total of 100,000 trials

Unconditional probability of detecting a twin peak: 0.00512 (0.51%)

99% CI: [0.00455, 0.0057]

The Clopper-Pearson method uses the Cumulative Distribution Function (CDF) of the binomial distribution to get precise bounds and is especially useful when dealing with a small number of successes. We use consistent units for the spectra throughout this analysis.

Simulation-2: In this simulation (instead of using 1/f decay), we start with the smoothed channel 5’s spectrum (Savitzky-Golay filter), which already has a dominant peak at frequency f₁=22.83Hz, and generate 40,000 surrogate trials. Similar to simulation_1, we add noise to this spectrum in each trial. We note that the prominence criterion of a peak is not reliant on the height of the peak; instead, it only considers how well the peak stands out compared with its immediate surroundings. Since the existing peak f₁ in the spectrum is quite robust (dominant) we want to ensure that the random noise peak f₂ we detect in the trials not only has prominence ≥ 0.04 but is also robust (high enough compared with f₁) to exhibit the beating phenomena on which our hypothesis is based, i.e., it (in conjunction with f₁) can produce a detectable and observable tremor frequency of |f₁−f₂| as seen in the video of the patient. Therefore, on top of the conditions for detecting peaks and twin peaks used in simulation-1, we add the following condition on the random peak’s height: 0.6f_1−height < f_2−height < 1.4f_1−height. This condition is satisfied by all the LFP twin peaks detected by our algorithm in this paper. Out of a total of 40,000 surrogate trials, we extract the following statistics:

743 successful trials out of a total of 40,000 trials.

Probability of detecting |f₁−f₂|≤[4.68, 5.76]Hz, i.e., a twin peak: 0.0185 (1.8%)

99% CI: [0.01688, 0.02038]

Based on the results of these simulations, we conclude that the dominant twin peak structure observed in our LFP data is unlikely to arise by chance due to random noise only.

In this paper, we provide strong support for our hypothesis comprising Equations 5a–c. Specifically, for the first time in the literature, we connect features of the beta band spectrum with the defining features associated with the tremor dynamics measured from the patient’s video, such as tremor frequency and tremor fluctuation time. This connection is between two spectra that lie in different frequency bands; beta-band STN LFPs lie in the 13-31 Hz range while the physical tremors of PD patients lie in the 4–8 Hz (Bain, 2002) range. Furthermore, beyond the hypothesis itself, we also use communication theory to theoretically explain this connection between the two spectra by asserting that the tremor frequency oscillations are due to the beating of two dominant frequencies in the beta band of the STN LFPs.

To summarize, the paper concludes that properties of a PD patient’s tremor such as frequency and tremor duration can be obtained from the patient’s subthalamic nucleus beta-band spectrum. Furthermore, we conclude that such a dominant twin peak structure is unlikely to arise by chance only.

The work in this paper sheds light on new possibilities for interpreting and utilizing DBS recordings of PD patients. Furthermore, the identification of beta-band dominant twin peaks, specifically in the lateral channels of directional leads, may suggest new directions in electrode design and application. This specificity may pave the way for enhanced precision in signal reading, possibly leading to more effective treatment delivery. The dominant Twin Peaks may also have implications for research in closed-loop DBS (Feldmann et al., 2022; Koeglsperger et al., 2021; Plate et al., 2021; Vissani et al., 2020) for PD patients. The utilization of Google’s MediaPipe AI (see text footnote 1) framework introduces a simple approach to measuring tremor characteristics without the need for complex multiple-camera motion capture systems.

While the match between the predictions of our hypothesis and actual tremor measurements demonstrated in this paper is remarkable, we acknowledge the limitations of a single-patient study, since it cannot address questions, to be addressed in future research, like: (i) Are dominant Twin Peaks detected by using the methods presented in this paper found in the STN DBS recordings of all PD patients exhibiting tremors? (ii) Is our hypothesis comprising Equations 5a–c valid for all PD patients with tremors? (iii) Can these results be extended to Essential Tremor (ET) patients? Therefore, after receiving IRB (Institutional Review Board) approval, we are preparing to apply our research methodology to a larger cohort of patients.