Real-time phase and amplitude estimation of neurophysiological signals exploiting a non-resonant oscillator

A recent advancement in the field of neuromodulation is to adapt stimulation parameters according to pre-specified biomarkers tracked in real-time. These markers comprise short and transient signal features, such as bursts of elevated band power. To capture these features, instantaneous measures of phase and/or amplitude are employed, which inform stimulation adjustment with high temporal specificity. For adaptive neuromodulation it is therefore necessary to precisely estimate a signal's phase and amplitude with minimum delay and in a causal way, i.e. without depending on future parts of the signal. Here we demonstrate a method that utilizes oscillation theory to estimate phase and amplitude in real-time and compare it to a recently proposed causal modification of the Hilbert transform. By simulating real-time processing of human LFP data, we show that our approach almost perfectly tracks offline phase and amplitude with minimum delay and is computationally highly efficient.


Introduction
Adaptive neuromodulation is an evolving concept in the field of deep brain stimulation (DBS), which encompasses techniques that aim at delivering electrical current in response to the identification of prespecified biomarkers (Krauss et al., 2021). These markers can be recovered from a diverse range of signals, such as peripheral, surface electroencephalogram or deep brain recordings (Bouthour et al., 2019). However, common to these modalities is the mostly transient nature of derived biomarkers. For example, short-lived bursts of elevated amplitude in the beta frequency range have been established as a surrogate marker for impaired motor performance in patients with Parkinson's disease (PD) (Lofredi et al., 2019). In first closed-loop DBS studies, instantaneous measures of beta amplitude were directly used as a feedback signal to trigger stimulation when beta estimates crossed a predetermined threshold Tinkhauser et al., 2017). An alternative approach is to use phase-specific stimulation to effectively suppress high-amplitude signal episodes (Cagnan et al., 2013). Accordingly, phase-locked adaptive neuromodulation based on peripheral tremor data (Schreglmann et al., 2021) or electrocorticography (ECoG) (McNamara et al., 2020) has recently been demonstrated. For both approaches it is necessary to deliver stimulation in close temporal proximity to biomarker detection to render adaptive neuromodulation strategies successful, which calls for reliable and fast causal estimation of phase and amplitude of neurophysiological signals.
Here we present and test a causal method for phase and amplitude computation in real-time. This technique exploits the computer model of a linear oscillator far from resonance and recomputes the phase and amplitude of this virtual oscillator into those of the analyzed signal (for mathematical details and proof of the principle, see Rosenblum et al. (2021)). We compare our approach against the offline, non-causal Hilbert Transform and its causal modification, endpoint-corrected Hilbert transform (ecHT) (Schreglmann et al., 2021) by simulating real-time processing of human local field potential (LFP) data. We demonstrate that our approach is more precise and much faster than ecHT.

Phase and amplitude measurement via a non-resonant oscillator
Our virtual "measurement device" exploits the property of a linear driven oscillator. If the input's frequency is much less than the resonance frequency of the oscillator, then the phase of the input approximately equals that of the oscillator, and the amplitudes differ by a known constant factor. Thus, practically, the problem boils down to solving numerically the equation of the oscillator driven by the signal of interest, ẍ + αẋ + ω 2 x = s(t k ), where α, ω are the oscillator's damping coefficient and frequency, respectively, and s(t k ) = s(k/f s ) is the input signal, sampled with the rate f s ; k = 1,2, … is the index. Next, for each time point and phase θ(t k ) = arctan(− αν/ω 2 − ν 2 ) and recompute them into input's amplitude a(t k ) and phase ϕ(t k ). Here ν is the frequency of s. The key idea is that for properly chosen oscillator parameters, the dependence on ν is weak so that ν shall be known only approximately. So, for ω ≫ ν and sufficiently small α we have ϕ(t k ) ≈ θ(t k ), while for larger α we compute factor. The described procedure is implemented by several lines of code, available in the supplementary materials.

Patients and recordings
7 patients (median age 55, range 44-65 years; 3 women) with PD who underwent DBS surgery targeting the subthalamic nucleus were included in this study. All patients provided written informed consent. The study was approved by the local ethics committee at Charité Universitätsmedizin Berlin (study number EA2/129/17) and followed the standards of the Declaration of Helsinki. Either the 4-contact electrode model 3389 (Medtronic, Minneapolis, USA) or 8-contact directional leads (Boston Scientific, Marlborough, USA) were implanted. LFP recordings were performed by a standardized procedure (Feldmann et al., 2021) with patients being assessed after at least 12 h withdrawal of dopaminergic medication. Patients were quietly seated and instructed to avoid any movement. Data were acquired from 8 (Medtronic) and 16 (Boston Scientific) electrode contacts, respectively, with a SAGA 64 amplifier (TMSi, Oldenzaal, The Netherlands) at 4000 Hz or 4096 Hz sampling rate and online referenced to the upper-or lowermost contact of one hemisphere.

Preprocessing
All analyses were carried out using custom-made scripts in Matlab 2020b (The MathWorks, Natick, USA). To allow for equivalent electrode layouts, LFPs from the 3-contact directional levels of the 8-contact directional electrodes were summed up to simulate a 4-contact layout, as previously reported (Lofredi et al., 2019). LFP data was then rereferenced bipolarly between adjacent contacts for each hemisphere. From each dataset, a continuous segment of 2 min duration was selected. Power spectral density estimates were obtained using Welch's method. A clear peak in the beta frequency range was identified in 13 hemispheres, while data from 1 hemisphere was excluded due to absent beta band elevation. The channel with the most prominent beta peak per hemisphere was selected. Before further analyses, LFP data were downsampled to a common sampling rate of 1 kHz and high-pass filtered with a cut-off frequency of 1 Hz.

Phase and amplitude estimation
For each dataset, phase and amplitude were computed with a nonrealizable "gold standard" technique (i.e. offline Hilbert transform (oHT)), with ecHT and with non-resonant oscillator (NRO), with the latter two being real-time capable. For oHT, before the Hilbert transform, LFP data was causally filtered using a 281 points FIR bandpass filter with the cut-off frequencies at ±3 Hz around the identified beta peak. Phase and amplitude estimates were then obtained by computing the instantaneous magnitude and phase angle of the analytic signal, respectively. For the other two methods, real-time processing was simulated by looping over samples of each dataset, using only the present and past samples. At each iteration, the signal was causally filtered as described above, followed by phase and amplitude estimation. For ecHT, the bandpass of the intrinsic Butterworth filter was set to ± half the frequency of the identified beta peak. Filter coefficients were computed outside the real-time loop. The analytic signal, as derived by ecHT, was computed at each iteration over a running window of 256 samples. Real-time phase and amplitude were recovered from the last sample of the running window. For NRO, ν was set to the frequency of the beta peak, ω was 5ν and α was 10 and 80 for phase and amplitude estimation, respectively. Upon each iteration, the most recent sample was fed to the algorithm. Processing time was measured by taking timestamps prior to and after the real-time simulation and dividing through the number of iterations.

Statistics
The first and last 5 s of each phase and amplitude trace were removed to avoid edge artifacts. Correlation of phase and amplitude estimates from both real-time techniques to oHT was computed by means of Pearson correlation coefficients for each dataset. To account for the discontinuous nature of phase, the cosine of phase estimates was computed beforehand. Paired sample t-tests were carried out on the Fisher z-transformed correlation coefficients for both phase and amplitude to assess performance differences between ecHT and the NRO. To investigate frequency-dependent variations in estimation performance, the Fisher z-transformed R-values were also correlated to the beta peak frequencies by means of Pearson correlation. The delay between both real-time techniques and oHT was determined by means of crosscorrelation. To test the possibility that a filter bandwidth other than suggested by the authors (Schreglmann et al., 2021) may improve ecHT performance, real-time phase and amplitude was recomputed with ecHT for dataset 6 employing different relative bandwidths. The set of phase and amplitude estimates obtained by this procedure were then correlated with oHT as laid out above.

Results
Median peak frequency within the beta range was 18 Hz (range 14.5 Hz -29 Hz, Fig. 1 A). Exemplary sections of phase and amplitude estimates are depicted in Fig. 1 B & C. Median correlation coefficients with oHT were 0.71 (ecHT, range 0.63-0.85) and 0.99 (NRO, 0.99-1) for phase and 0.66 (ecHT, 0.6-0.88) and 0.99 (NRO, 0.99-1) for amplitude, respectively ( Fig. 1 D & E). Paired sample t-tests on the Fisher z-transformed correlation coefficients revealed significant differences between ecHT and NRO for both phase (t(12) = − 132.2, p < 0.001) and amplitude (t(12) = − 73.1, p < 0.001). Both real-time techniques showed strong frequency-dependent performance variations for both phase and amplitude estimation with beta peak frequency positively correlating with estimation quality (all R > 0.83, all p < 0.001). Median delay between oHT and the real-time methods were 57 ms (ecHT, 35 ms -70 ms) and 0 ms (NRO, 0 ms -0 ms) for phase and 54 ms (ecHT, 32 ms -67 ms) and 0 ms (NRO, 0 ms -1 ms) for amplitude, respectively. The relative bandwidth sweep for dataset 6 yielded maximum correlation coefficients of 0.87 and 0.93 at relative bandwidths of 1.25 and 1.5 for phase and amplitude, respectively. Median computation time for both amplitude and phase was 0.02 ms (ecHT) and 0.0006 ms (NRO) on an iMac (Intel Quad-Core i5 processor, 3.8 GHz).

Discussion
In this report we present a real-time capable technique to estimate amplitude and phase with minimum delay. We compared this approach to a recently suggested method (ecHT) by correlating the output from both approaches to a non-realizable "gold standard" (oHT). Additionally, we determined computation times of both methods.
Our method significantly outperformed ecHT in estimating both phase and amplitude of high-frequency signals. For phase, NRO reliably tracked oHT, while ecHT failed to provide reasonable results for some datasets. Furthermore, while changing filter parameters of ecHT improved estimation performance to some extent, NRO still yielded superior results. With a median delay of 57 ms (ecHT) versus 0 ms (NRO) and a cycle length of 55 ms (translating to a frequency of 18 Hz), phase estimates derived from ecHT lagged behind oHT more than a full oscillatory cycle.
Performance differences for amplitude estimation were similarly pronounced. Note that ecHT generates a smoother amplitude trace, which may be beneficial for threshold-based adaptive neuromodulation in that stimulation amplitude fluctuates to a lesser extent. Indeed, pilot studies of adaptive DBS in PD utilized strongly smoothed amplitude time series Tinkhauser et al., 2017). However, which amount of smoothing proves to be most effective clinically is an open question.
Correlating beta peak frequencies with estimation performance revealed that for both techniques reduced accuracy was associated with lower frequencies. However, while these frequency-dependent differences were negligible for NRO, considerable performance drops were observed for ecHT. This may be of particular relevance for beta-driven closed-loop DBS for PD as elevated band-power specifically in the low beta frequency range is indicative for motor symptom severity (Neumann et al., 2016).
Importantly, increased beta power in PD is not a stationary phenomenon, but rather occurs in phasic bursts which accompany motor deterioration (Lofredi et al., 2019). Thus, adaptive DBS that addresses dynamic fluctuations in beta power as a biomarker may prove beneficial by disrupting a putative pathophysiological mechanism (Little and Brown, 2020). Recent closed-loop DBS implementations have therefore employed instantaneous measures to capture transient states of hypersynchronization. In these studies, high-frequency stimulation was either triggered by exaggerated amplitude itself Tinkhauser et al., 2017), or stimulation pulses were delivered at specific oscillatory phases aiming for stronger desynchronizing effects (McNamara et al., 2020). Here, NRO was shown to provide highly accurate measures of instantaneous phase and amplitude in real-time, which may therefore substantially improve the timely discharge of stimulation within an adaptive neuromodulation device.
Common to both techniques is the issue that any upstream causal filter induces additional delay per se. We have excluded this aspect from our analysis by causally filtering the signals prior to applying both offline and real-time techniques. Thus, our analyses are only informative in regard to any further delay that is introduced by the algorithm used for phase or amplitude estimation. However, more elaborate techniques are available to reduce filter delay (Smetanin et al., 2020), which may be used in conjunction with our method. Similarly, signal quality might Fig. 1. A: Power spectra of all datasets. Columns indicate patients and rows indicate hemispheres (first row = right hemisphere, second row = left hemisphere). Selected peak frequencies are marked by vertical bars in dark grey and filter bandwidths are indicated by areas around each peak frequency in light grey. B, C: Exemplary sections from dataset 1 of estimated phase (B) and amplitude (C) for all techniques employed. D, E: Correlation coefficients between oHT and ecHT (red) and between oHT and NRO (blue), for phase (D) and amplitude (E) estimates for each dataset.
impact estimation performance. To this end, we high-pass filtered each dataset prior to conducting phase and amplitude estimation. While this procedure is non-causal per se, we have adopted it in our simulations to rule out performance drops due to technical artifacts unrelated to any of the techniques employed.
Computation time for our algorithm was around 30 times faster than ecHT. However, absolute processing times were still far below the sampling frequency employed here for both techniques. Thus, both methods are real-time capable not only theoretically, but also in practice. However, depending on downstream computations performed after phase or amplitude estimation within an adaptive neuromodulation device, computational advantages might still hold beneficial . Additionally, in practical applications further sources of computational delay have to be considered as well (e.g. communication time between an external processor and the neurostimulator), calling for fully integrated, embedded systems (Opri et al., 2020).
To summarize, we have described a new method for real-time estimation of phase and amplitude that induces minimum delay and is computationally efficient. This method may be adopted in adaptive neuromodulation devices to quickly detect biomarkers based on amplitude or to track oscillatory phase, supporting the timely delivery of brain stimulation.

Declaration of Competing Interest
AAK declares that she is on the advisory board of Boston Scientific and Medtronic, and has received honoraria from Boston Scientific and Medtronic.