Removal of Electrocardiogram Artifacts From Local Field Potentials Recorded by Sensing-Enabled Neurostimulator

Sensing-enabled neurostimulators are an advanced technology for chronic observation of brain activities, and show great potential for closed-loop neuromodulation and as implantable brain-computer interfaces. However, local field potentials (LFPs) recorded by sensing-enabled neurostimulators can be contaminated by electrocardiogram (ECG) signals due to complex recording conditions and limited common-mode-rejection-ratio (CMRR). In this study, we propose a solution for removing such ECG artifacts from local field potentials (LFPs) recorded by a sensing-enabled neurostimulator. A synchronized monopolar channel was added as an ECG reference, and two pre-existing methods, i.e., template subtraction and adaptive filtering, were then applied. ECG artifacts were successfully removed and the performance of the method was insensitive to residual stimulation artifacts. This approach to removal of ECG artifacts broadens the range of applications of sensing-enabled neurostimulators.


INTRODUCTION
Recently, sensing-enabled neurostimulation has emerged as a technology for long-term observation of brain activities and has paved the way for development of closed-loop neuromodulators and implantable brain-computer interfaces (Stanslaski et al., 2012;Qian et al., 2014;Martini et al., 2020;Ramirez-Zamora et al., 2020). Using sensing-enabled neurostimulators, a seminal series of studies made enormous advances in the mechanisms of deep brain stimulation (DBS) (Trager et al., 2016;Neumann et al., 2017), the effects of closed-loop neuromodulation (Meidahl et al., 2017;Velisar et al., 2019), and the feasibility of implantable brain-computer interfaces (Vansteensel et al., 2016;Golshan et al., 2020). By integrating recordings of LFPs, sensing-enabled neurostimulators have opened a real window into chronic brain activities (Shen, 2014). However, it is difficult to faithfully record LFPs due to electrocardiogram (ECG) and stimulation artifacts (Sorkhabi et al., 2020). Previous studies using sensing-enabled neurostimulators reported considerable data loss due to ECG contamination (Quinn et al., 2015;Trager et al., 2016;Anidi et al., 2018;Hell et al., 2018;Swann et al., 2018). The research of Quinn et al. (2015) showed that six of the sixteen recordings were excluded because of ECG artifacts. In some cases, ECG artifacts could develop in the LFPs recorded during longitudinal followups (Trager et al., 2016). The researchers had to choose those LFP channels free of ECG contamination (Swann et al., 2018), which limited the precision of the targeted recording positions.
ECG artifacts can be attributed to an inadequate commonmode-rejection-ratio (CMRR) of the recording module in a sensing-enabled neurostimulator. For LFP recording, signals were usually differentiated between pairs of contacts on the electrode. The ECG spikes are regarded as commonmode signals that can be rejected by differentiating. To eliminate ECG artifacts, previous studies have shown that the CMRR has to be greater than 60 dB (Sorkhabi et al., 2020). However, achieving a high and stable CMRR is very challenging for implantable devices because power consumption and size of the neurostimulators are very limited. Furthermore, after implantation, a slight leakage of fluid into the neurostimulator can break the symmetry of the differentiated channels, which alters the CMRR (Quinn et al., 2015).
Although there are ample algorithms for removal of ECG artifacts from electrophysiological signals, few can be applied to LFP recordings. Template subtraction calculates the average waveform of ECG artifacts and subtracts it from each spike in contaminated recordings (Zhou et al., 2007;Marker and Maluf, 2014). A previous study used a raw ECG signal and template subtraction to remove the ECG artifacts from LFPs recorded by sensing-enabled neurostimulators (Canessa et al., 2016). Adaptive filtering is another popular method for removing ECG artifacts (Lu et al., 2009;Sweeney et al., 2012). It estimates the noise component and subtracts it from the original recording. Both methods usually require a synchronized ECG waveform for reference. Independent component analysis (ICA) is a widely used blind-source separation method that dispenses with the need for an ECG reference (Mak et al., 2010;Sweeney et al., 2012). However, the recordings of sensing-enabled neurostimulators are usually made without an ECG reference, and the number of LFP channels is too limited (for example, two channels) to perform ICA. Therefore, these methods are unsuitable for removal of ECG artifacts from most LFP recordings. Although there are some blind-source separation methods for single-channel recordings (Sweeney et al., 2012), their robustness is inconclusive due to residual DBS artifacts.
In this study, we propose a method for the removal of ECG artifacts from LFPs. A sensing-enabled neurostimulator was used to record simulated LFPs contaminated with ECG artifacts in saline solution. We modified the recording montages to add a synchronized ECG reference channel. Using the reference channel, we used template subtraction and adaptive filtering to remove ECG artifacts from simulated LFPs. We evaluated the performance of the proposed method of ECG artifact removal and explored the effects of DBS artifacts. The results show that using modified recording montages allows ECG artifacts to be successfully removed, thus revealing the LFP signals.

Sensing-Enabled Neurostimulator
The LFP recording instrumentation is illustrated in Figure 1. A sensing-enabled neurostimulator (G102RS, Beijing PINS Medical Co., Ltd.) which can be fully implanted for DBS therapy was used. Before sampling, signals were first filtered by a builtin 0.3-250 Hz band-pass filter. Then the signals were sampled by the neurostimulator and wirelessly transmitted to the external PC platform through radio frequency communication. The transmission rate was 250 kbps and the delay was less than 10 ms. The wireless communication distance was approximately 2 m. A rechargeable battery was incorporated into the neurostimulator to ensure long-term recording.

Modification of Recording Montages
For LFP recording, signals were usually differentiated between pairs of contacts in the electrode, i.e., bipolar recording. For example, in Figure 2, the signal was differentiated between contact 1 and contact 3 in the DBS electrode after the passive filter FIGURE 1 | Illustration of the sensing-enabled neurostimulator system. FIGURE 2 | Schematic of the recording module in the sensing-enabled neurostimulator. LFPs were recorded in the bipolar recording channel (between contacts 1 and 3). The ECG signals were recorded in the monopolar recording channel (between the titanium case and contact 4). The capacitors and resistors in each chain formed passive band-pass filters (0.3-250Hz) of the analog-front-end circuits. Monopolar stimulation was delivered between the titanium case and contact 2. network. Monopolar stimulation was synchronously delivered between the titanium case of the neurostimulator (anode, implanted in the chest) and the therapeutic contact 2 (cathode). Potentials generated by the stimulation and the ECG source were approximately the same at contacts 1 and 3. Thus, both of the ECG potentials and stimulation pulses could be regarded as common-mode signals in the bipolar recording chain.
To achieve efficient ECG artifact removal, we developed recording montages by adding a channel for the ECG reference in the neurostimulator. As shown in Figure 2 (red lines), a monopolar recording channel was differentiated between the titanium case and another contact in the DBS electrode. In this monopolar channel, ECG signals were synchronously recorded with the LFPs in the bipolar recording channel. Although the monopolar recording montage may yield larger DBS artifacts, recording ECG signals is possible as long as the analog-to-digital converter is not saturated. In the following sections, we describe tests of the feasibility of synchronized bipolar and monopolar recordings during stimulation and evaluate the benefits of the modified recording montages for removal of ECG artifacts.

Simulated Recording Setup
As shown in Figure 3, to simulate signal recordings after implantation, the neurostimulator's titanium case and the DBS electrodes were fully immersed in an ASTM phantom (ASTM, 2011) filled with saline solution at room temperature. The DBS electrodes (model L301, Beijing PINS Medical Co., Ltd.) have four platinum-iridium cylindrical contacts. The contacts were 1.3 mm in diameter, 1.5 mm in length, and spaced 0.5 mm apart. Two Ag/Cl disc-electrodes beside the DBS FIGURE 3 | Illustration of the simulated recording setup. The titanium case and the DBS electrodes were fully immersed in the phantom that was filled with saline solution.Phantom size is expressed in centimeters. A 23Hz sinusoidal signal and a standard ECG waveform were delivered by two pairs of Ag/Cl disc-electrodes. The recordings were wirelessly transmitted toa computer. electrodes delivered a 23 Hz sinusoidal signal that simulated an LFP, which was chosen for the following reasons: (1) the beta band (13-35 Hz) of LFPs recorded by sensingenabled neurostimulators is the most relevant band for revealing the mechanisms of DBS and for developing a closed-loop stimulation strategy (Trager et al., 2016;Meidahl et al., 2017;Neumann et al., 2017;Velisar et al., 2019), and (2) based on previous clinical experience, ECGs mainly contaminate the low-frequency band (below 100 Hz) in these specific LFPs, especially in the beta band. For the relevant signals in the high-frequency band (above 100 Hz), a high-pass filter could remove the ECG artifacts (Redfern et al., 1993;Zhou and Kuiken, 2006;Marker and Maluf, 2014). The method proposed in the current study aims to extract LFPs in the ECGcontaminated frequency band, i.e., the low-frequency band. Thus, a 23 Hz tone should be sufficient for validation of the method. Two additional Ag/Cl disc-electrodes beside the titanium case delivered a standard ECG waveform generated from a digital function generator. The period of the ECG signal was set to 750 ms.

Stimulation and Recording
Electrode contacts 1 and 3 were configured as the bipolar montage for LFP recording. The sensing-enabled neurostimulator case and electrode contact 2 were configured as the monopolar montage for ECG reference recording. The bipolar and monopolar recordings were synchronized by the sensing-enabled neurostimulator. The sampling rate was 1,000 Hz. Monopolar stimulation was delivered between the case and contact 2 in the DBS electrodes. In a real recording scenario after the implantation of the sensing-enabled neurostimulators, no power line noise could be sensed. However, in the saline simulation, strong power line noises could be sensed by the neurostimulator. To avoid cross-interference (differential harmonics) between stimulation and the power frequency (50 Hz), the stimulation frequency was set to 150 Hz. The stimulation pulse width was 60 µs. To study the effects of DBS artifacts, the stimulation amplitude was set to off, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 V. A segment of clean LFP was recorded after the ECG source and stimulation were both turned off. Under each condition, data was recorded for at least 300 s.

Data Preprocessing
All recordings were first filtered by a 1 Hz high-pass filter (digital Butterworth) to suppress the baseline fluctuations and then filtered by a 470 Hz low-pass filter (digital Butterworth) to exclude higher aliasing frequencies. For recordings of a clean LFP, the power frequency of 50 Hz and its harmonics at 100 Hz were filtered by an adaptive trap filter (Keshtkaran and Yang, 2014). For the recordings of ECG-contaminated LFP, DBS artifacts combined with the power frequency inferences were filtered out. To retain the information in the frequency domain as much as possible, we used a custom-designed adaptive trap filter to remove the DBS artifacts. For recordings of the ECG reference signal, data was filtered by the adaptive filter used for the power frequency and then filtered by a 100 Hz low-pass filter (3rd-order digital Butterworth).

ECG Artifact Removal
In this study, a monopolar channel was added for ECG recording. We took advantage of two pre-existing methods, i.e., template subtraction and adaptive filtering, to evaluate the effectiveness of ECG artifact removal. MATLAB (version 2016a) was used for all the calculations.

Template Subtraction
The template subtraction method extracts an ECG artifact template from the contaminated LFP recordings and subtracts the template in each ECG artifacts spike. The method mainly involved three steps: (1) ECG spike detection, (2) ECG artifact template extraction, and (3) ECG artifact template subtraction. For ECG spike detection, we applied an accurate QRS complex detector proposed by Pan and Tompkins (1985) and Hamilton and Tompkins (1986) to the ECG reference channel. The detected R-waves were used to align the ECG artifacts in the contaminated LFPs. A template was extracted by averaging the ECG artifacts and filtered by a 100 Hz low-pass filter (3rd-order digital Butterworth). The filtered template was then subtracted from the ECG-contaminated LFPs.

Adaptive Filtering
The general methodology of adaptive filtering is to estimate the ECG artifacts from the contaminated LFPs by minimizing the errors between the output and the ECG reference. In this study, we took advantage of the built-in frequency domain adaptive filter toolbox (dsp.FrequencyDomainAdaptiveFilter) in MATLAB to design the adaptive filter. The length of the coefficient vector of the adaptive filter was 64 and the block length of the coefficients updates was 1,000. The frequency-domain adaptive filter iterated the coefficients vector after transformation of the inputs to the frequency domain using discrete Fourier transforms. For detailed algorithms, please see Shynk (1992). The input was the ECG reference and the desired signal was the contaminated LFP. The outputs were the estimated ECG artifacts. The estimated ECG artifacts were then subtracted from the contaminated LFP recordings.

Performance Evaluation
Power spectral density (PSD) was estimated by the Welch's method. The frequency resolution of the PSDs was 1 Hz and the overlap of the Hamming window was 50%. Logarithmic PSDs were plotted for comparison between the clean LFP, the contaminated LFPs, and the cleaned LFPs. To quantitatively evaluate artifact removal, the signal-noise-ratio (SNR) and the root mean square logarithmic error (RMSLE) of the PSDs were calculated.
The SNR was defined as the ratio of the average power in the signal and noise bands. The SNR was calculated as: where Power signal represents the average power of the signal in the band of 22-23 Hz, Power noise represents the average power of the ECG contamination in the band of 1-100 Hz, excluding 22-23 Hz. A larger SNR indicates better ECG artifact removal.
Because the spectrum was not uniformly distributed across the band, we used the RMSLE to evaluate the difference between the PSDs of cleaned LFPs and the simulated clean LFP. By logarithmic transformation, RMSLE evaluates the differences between small values (such as the power of the noise) similarly to differences between big values (such as the power of the signal). The RMSLE was calculated as: where psd sin (f) represents the logarithmic PSD of the clean LFP, psd(f) represents the logarithmic PSD of the contaminated or cleaned LFPs. f 1 and f 2 indexes the relevant band for analysis, which is 1-100 Hz (i.e., N = 100). A smaller RMSLE indicates better performance of ECG artifact removal.

Simulated LFP Recordings
In this study, we synchronously recorded bipolar and monopolar channels for ECG artifact removal. Examples of the recordings are shown in Figure 4A. The top axis in Figure 4A is the time series of bipolar recordings. The dark line is the clean LFP and the gray line is the contaminated LFP. The residual DBS (2.5 V, 150 Hz, 60 µs) artifacts in the contaminated LFP are more than two orders of magnitude greater than in the clean LFP. The ECG artifacts could hardly be distinguished in the time series. The bottom axis in Figure 4A is the time series of the monopolar recording synchronized with the contaminated LFP. Although the DBS artifacts were greater in the monopolar recording, a clear ECG waveform could be distinguished. Figure 4B shows the PSD of the clean LFP. Figure 4C shows the PSD of the contaminated LFP. The fundamental frequency of the DBS artifacts was 150 Hz and the harmonic frequencies were 300 and 450 Hz. Other aliasing frequencies superposed with the power frequency were integer multiples of 50 Hz.

ECG Recording Performance
To extract clear ECG spikes, monopolar recordings were filtered by a fixed 100 Hz low-pass filter (3rd-order digital Butterworth).
To study whether the ECG spikes in the monopolar recordings could be consistently restored in the presence of different DBS artifacts, we compared the amplitudes of the extracted spikes. The spikes were aligned according to the detected R-waves. Figure 5A shows the averaged spikes extracted under different stimulation amplitudes. Figure 5B shows the peak-to-peak amplitudes of each ECG spike and the peak-to-peak amplitudes of the DBS artifacts in the original monopolar recordings. Compared with the mean value in the DBS off state (5.57 ± 0.60 mV, mean ± SD), the amplitudes of the spikes decreased slightly in the DBS on state (4.54 ± 0.66 mV, p < 0.001, N = 266). Although the amplitudes of the DBS artifacts increased from 2.50 ± 1.8 mV (0.5 V DBS) to 9.65 ± 0.29 mV (3.0 V DBS), no significant changes were found  between the amplitudes of the ECG spikes (p > 0.064, N = 266). The results indicated that by using a simple fixed 100 Hz lowpass filter, ECG spikes could be extracted from the monopolar recordings under different stimulation amplitudes. The robust performance of spike extraction provided a foundation for ECG artifact removal.

Performance of ECG Artifact Removal
To validate the feasibility of using the monopolar channel as a reference for ECG artifact removal, we tested two preexisting methods, i.e., template subtraction and adaptive filtering. Performance of ECG artifact removal using each method was evaluated in the DBS off and on states. Figure 6 shows the results of template subtraction. In Figure 6A, the red line on the top axis is the ECG-contaminated LFP recorded in the DBS off state and the gray line is the cleaned LFP. The bottom axis shows the distinct R-waves in the synchronized monopolar channel. Using the monopolar channel as an ECG reference, the template subtraction method removed the ECG artifacts in the contaminated LFP. Figure 6B shows the PSDs of the ECG-contaminated LFP (red line) and cleaned LFP (gray line). The ECG artifacts mainly contaminated the band below 100 Hz. After removing the artifacts, the PSD of the cleaned LFP was overlaid with the PSD of the simulated clean LFP (blue line). In Figure 6C, an ECGcontaminated LFP recorded in the DBS on state (2.5 V, 150 Hz, 60 µs) and the corresponding cleaned LFP are plotted on the top axis. After removing the DBS artifacts, the monopolar recording reflected distinct R-waves (bottom axis). Using the monopolar channel as an ECG reference, the ECG artifacts were successfully removed. Figure 6D shows the PSDs of the contaminated (red line) and cleaned (gray line) LFPs in the DBS on state. The ECG-contaminated band below 100 Hz was restored in the cleaned LFP. We found a slightly elevated noise floor in the PSD of the cleaned LFP, which partially resulted from the ECG signal generator. In addition, some noise came from residual stimulation artifacts. Although the stimulation artifacts were largely removed, some harmonics and aliasing components remained, such as at 250, 300, 400, and 450 Hz. Spectrum estimation of these residual stimulation artifacts would unavoidably introduce leakage of the power to a broadband range. Figure 7 shows the results of adaptive filtering. Figures 7A,B show recordings in the DBS off state, and Figures 7C,D show recordings in the DBS on state (2.5 V, 150 Hz, 60 µs). The ECGcontaminated band below 100 Hz was restored in the cleaned LFPs. Therefore, both template subtraction and adaptive filtering could remove the ECG artifacts.
To quantitatively evaluate the performance of the method, SNRs and RMSLEs of the PSDs were calculated. Table 1 shows the SNR and RMSLE of the cleaned LFPs. Compared with the ECG-contaminated LFPs, the SNRs of the cleaned LFPs were greatly improved, with values close to the SNR of the clean LFP (20.87 dB). The RMSLEs also decreased, indicating that the differences in the PSDs between the cleaned LFP and the simulated clean LFP decreased after ECG artifact removal. Neither SNR nor RMSLE were significantly influenced by increasing the stimulation amplitude.

Performance Under Residual DBS Artifacts
In section "ECG Artifact Removal, " the DBS artifacts under different stimulation amplitudes were precisely filtered. However, due to complex recording conditions and various filtering methods, DBS artifacts could remain in the LFPs at varying levels. Therefore, it is crucially important to evaluate the performance in a background of residual DBS artifacts. Figure 8A shows the contaminated LFP in the stimulation on state (gray lines, 2.5 V, 150 Hz, 60 µs) without filtering DBS artifacts. The recordings were aligned according to the ECG spikes detected in the monopolar channel. By overlapping the recordings, the envelope of the ECG artifacts could be distinguished (dashed red lines). The white line in Figure 8A shows the average waveform of the overlapped recordings. Using this template, the ECG artifacts were subtracted from the contaminated LFP. Figure 8B shows the PSDs of the contaminated LFP (red line) and ECGcleaned LFP (gray line). The ECG-contaminated band was restored while the DBS artifacts remained. Figure 8C shows the ECG templates extracted from the contaminated LFPs. The performance of template extraction was stable across 0.5-3.0 V DBS without additional filtering of the DBS artifacts. Table 2 shows a quantitative evaluation of the ECG removal performance without filtering the DBS artifacts. To exclude the effects of increasing the DBS artifacts, SNR and RSMLE were calculated for the band below 40 Hz. Compared with The values in the brackets are the differences with the contaminated local field potentials (LFPs). # Thedeep brain stimulation (DBS) artifacts in the recordings were removed. * TS, template subtraction; AF, adaptive filtering.
the ECG-contaminated LFPs, both the SNR and RMSLE of the cleaned LFPs were improved. SNR of the cleaned LFPs was close to the value of the simulated clean LFP (18.70 dB). From 1.0 V DBS to 3.0 V DBS, the SNR of the cleaned LFPs decreased from 18.29 to 17.62 dB and the RSMLE increased from 0.98 to 1.9 dB, indicating an attenuation of performance. However, this attenuation of the performance was very slight compared with the increase on the DBS artifacts from 1.0 V DBS (211.01 µV) to 3.0 V DBS (516.18 µV). The results indicate that the ECG artifacts in the contaminated LFPs can be suppressed to a significant extent even in the presence of residual DBS artifacts.  Signal-noise-ratio (SNR) androot mean square logarithmic error (RMSLE) were calculated for the band below 40 Hz. The values in the brackets are the differences with the contaminated local field potentials(LFPs).

DISCUSSION
In this study, we proposed an approach to removing ECG artifacts from the signals recorded by sensing-enabled neurostimulators. A simultaneous monopolar montage was added and distinct R-waves could be recorded in the time series. Using the monopolar recording as an ECG reference, conventional ECG filtering methods, including template subtraction and adaptive filtering, could effectively remove ECG artifacts in the bipolar LFP recordings. The performance of this method was not significantly influenced by the magnitude of residual DBS artifacts. Combining the monopolar and bipolar recordings, the method could remove ECG artifacts without additional complex signal processing.

Using Simultaneous Monopolar Recording as an ECG Reference
Removing ECG artifacts is a common challenge in the field of electrophysiological signal processing. There have been many algorithms proposed in previous studies, such as template subtraction (Zhou et al., 2007;Marker and Maluf, 2014), adaptive filtering (Lu et al., 2009;Sweeney et al., 2012), ICA (Mak et al., 2010), and other blindsource separation methods (Sweeney et al., 2012). Most of the methods proved effective for multi-channel EEG and electromyography recording. Some methods were specifically designed for single-channel recording (Sweeney et al., 2012). However, most of the pre-existing methods cannot be directly applied to the signals recorded by implanted sensing-enabled neurostimulators. For sensing-enabled neurostimulators, the electrophysiological signals (LFPs) are usually independently recorded without ECG reference. Thus, the methods of template subtraction or adaptive filtering, both of which require reference signals, are not suitable. The number of synchronized channels of sensing-enabled neurostimulators is usually very limited. For recording in the DBS on state, the sensing-enabled DBS usually records only a few synchronized channels due to limited electrodes. Thus, the method of ICA which needs multi-channel recordings is not fully applicable to the LFPs recorded during DBS. Although there are some single-channel blind-source separation methods (Sweeney et al., 2012), the cost of the algorithms is usually higher and the robustness of the performance is inconclusive due to residual DBS artifacts. Lack of an ECG reference and recording channels makes it difficult to remove ECG artifacts in the LFPs recorded by sensing-enabled neurostimulators.
There are two kinds of recording montages in sensingenabled neurostimulators, i.e., the bipolar montage and the monopolar montage. The bipolar montage differentiates the potentials between pairs of contacts in the DBS electrode while the monopolar montage differentiates the potentials between one contact in the DBS electrode and the titanium case of the neurostimulator (Stanslaski et al., 2012;Qian et al., 2014Qian et al., , 2017. Previous studies reported that the monopolar montage is considerably affected by remote volume conductivity (Marmor et al., 2017). Thus, it is possible to record ECG references using a monopolar montage.
In this study, we demonstrated the feasibility of recording ECG reference signals using the monopolar montage of a sensingenabled neurostimulator. Our results showed that although the DBS artifacts were large, a clear R-waveform could be extracted in the monopolar channel. However, the peak-topeak amplitudes of ECG artifacts in bipolar recordings in vivo might be unstable due to respiration or motion. Previous literature reported that respiration or even relatively slow motion over time would induce baseline wander in the ECG signals (Satija et al., 2018;Chatterjee et al., 2020). The peak-to-peak amplitudes of ECG spikes vary by 15% due to baseline wander (Friesen et al., 1990;Satija et al., 2018). For rapid motion artifacts, the effect is usually over a short duration and the fluctuation could be 500 percent of the peak-to-peak amplitudes of ECG spikes (Friesen et al., 1990;Satija et al., 2018). In Figure 6C, the coefficient of variation of the ECG peak-topeak amplitudes in bipolar recordings was 13.07% (std/mean) and the value of the ECG spikes in monopolar recordings was 0.7%. To estimate the influence of these variations, we added amplitude-modulated white noise to the recordings. The coefficients of variation increased to 21.62% in bipolar recordings and 16.37% in monopolar recordings. Both variations in the bipolar and monopolar recordings were above the typical value reported in the literature. The methods of template subtraction and adaptive filtering were applied to remove ECG-artifacts in the noisy recordings. The results showed that the ECG artifacts could still largely be removed (see Supplementary Material for details).
The primary benefit of the current study is to make conventional ECG removal methods that require ECG reference signals applicable to LFPs recorded by sensingenabled neurostimulators. The method proposed in our study may help researchers to extract real LFPs to study chronic deep brain neural activities and DBS mechanisms. Most such research is conducted retrospectively. Where realtime analysis is needed, such as in closed-loop stimulation, a template could only be extracted in a short window. In a previous study using the template subtraction method for real-time ECG artifact removal, 30 ECG spikes were used to update the template (Abbaspour and Fallah, 2014). We tested time windows of 30, 20, and 10 s, respectively, and found that the template was valid in the 30 s window (see Supplementary Material for details). Because the interval of the simulated ECG spikes was 0.75 s, the results indicate that at least 40 ECG spikes would be necessary for extracting a valid template. In real-time scenarios, the templates could be accumulated and the performance could be improved by introducing overlap and padding techniques for each window. Alternatively, we could use the method of adaptive filtering for real-time analysis. Future studies should focus on development and refinement of detailed real-time algorithms based on specific real-time scenarios.

Effects of DBS Artifacts
The magnitudes of the stimulation pulses in the monopolar recording chain were asymmetric. The asymmetry would induce large DBS artifacts in the ECG reference channel. Consequently, different magnitudes of residual DBS artifacts might influence detection of ECG spikes. Previous studies have shown that as long as the amplifiers are not saturated, the real signals in an asymmetric montage can still be restored (Kent and Grill, 2012;Qian et al., 2017). In our study, we compared the amplitudes of the ECG spikes under different DBS amplitudes. The results showed that although the DBS artifacts increased greatly from 0.5 V DBS to 3.0 V DBS, no significant difference was found between the amplitudes of extracted ECG spikes. Robust R-waves could be extracted from the monopolar recordings. This result provided the foundation for the removal of ECG artifacts.
In the bipolar channel of the sensing-enabled neurostimulator, the amplitudes of DBS artifacts could be more than two orders of magnitude greater than the ECG artifacts and LFP signals. In our study, we precisely filtered the DBS artifacts to restore the time series. However, residual DBS artifacts were usually unavoidable in LFP recordings due to complex recording conditions and different filtering methods. To evaluate the effects of residual DBS artifacts on the removal of ECG artifacts, we studied the performance of the method without extra filtering of the DBS artifacts. The results showed that even in the presence of large DBS artifacts, ECG artifacts were suppressed significantly, indicating that the method was not sensitive to residual DBS artifacts. This finding emphasizes the robustness of the method and limits restrictions on its algorithms.

Limitations and Future Work
Removal of ECG artifacts in our study was only tested using the conventional fixed-template subtraction method and a simple adaptive filter. The performance could be improved by optimizing the extraction template and design of the adaptive filter. The current method was only tested on simulated signals recorded in saline solution. Although we have tested and verified the methods using physiological ECG levels, the performance of the method could still be influenced by chronic changes in the recording impedance, leakage of fluid, and different in-vivo electromagnetic environments. The performance needs to be confirmed for real in-vivo recordings in future studies. In addition, we should estimate whether the algorithms can be integrated into an embedded system. The goal of the current study was to provide an ECG reference for use with current methods and thus to reduce the complexity of the related algorithms. Based on our survey of the computation efficiency and storage resource of the current sensingenabled neurostimulators, both of the approaches proposed in this study could be completed within 400-500 ms. This is sufficient for real-time analysis. Thus, this method holds potential for development of algorithms for real-time removal of ECG artifacts in an embedded system. However, there is still plenty of work that needs to be done in future embedded development, especially for addressing the tradeoff between power consumption and analytical performance in specific applications.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation, to any qualified researcher.