Automatic removal of high-amplitude stimulus artefact from neuronal signal recorded in the subthalamic nucleus

Journal of Neuroscience Methods 198 (2011) 135–146

Contents lists available at ScienceDirect

Journal of Neuroscience Methods journal homepage: www.elsevier.com/locate/jneumeth

Automatic removal of high-amplitude stimulus artefact from neuronal signal recorded in the subthalamic nucleus Tarik Al-ani a,b,c,∗ , Fanny Cazettes c , Stéphane Palfi d , Jean-Pascal Lefaucheur e,f a

LISV-UVSQ, 10-12 Av. de l’Europe, 78140 Vélizy, France Dept. Informatique, ESIEE-Paris, Université Paris Est, Cité Descartes-BP 99, 93162 Noisy-Le-Grand, France ISBS-ESIEE-Paris, Cité Descartes-BP 99, 93162 Noisy-Le-Grand, France d Unité de Neurochirurgie Fonctionnelle, Hôpital Henri Mondor, APHP, 51 av. de Lattre de Tassigny, 94010 Créteil, France e Service de Physiologie, Hôpital Henri Mondor, APHP, 51 av. de Lattre de Tassigny, 94010 Créteil, France f EA 4391, Faculté de Médecine, Université Paris Est Créteil, 8 av. du Général Sarrail, 94010 Créteil, France b c

a r t i c l e

i n f o

Article history: Received 28 November 2010 Received in revised form 11 March 2011 Accepted 26 March 2011 Keywords: Stimulus artefact Deep brain stimulation Motor cortex stimulation Ensemble empirical mode decomposition Parkinson’s disease

a b s t r a c t Over the last few years, deep brain stimulation (DBS) with targets such as the subthalamic nucleus or the pallidum were found to be beneficial in the treatment of Parkinson’s disease and dystonia. The investigation of the mechanisms of action of DBS by recording concomitant neural activities in basal ganglia is hampered by the large stimulus artefacts (SA). Approaches to remove the SA with conventional filters, or other conventional digital methods, are not always effective due to the significant overlap between the spectral contents of the neuronal signal and the SA. Thus, such approaches may produce a significant residual SA or alter the neuronal signal dynamics by removing its frequency contents. In this work, we propose a method based on an on-line SA template extraction and on the Ensemble empirical mode decomposition (EEMD) to automatically detect and remove the dynamics of the SA without altering the embedded dynamics of the neuronal signal during stimulation. The results, based on real signals recorded in the subthalamic nucleus during Motor cortex stimulation (MCS) experiments, show that this technique, which may be applied on-line, effectively identifies, separates and removes the SA, and uncovers neuronal potentials superimposed on the artefact. © 2011 Elsevier B.V. All rights reserved.

1. Introduction Deep brain stimulation (HF-DBS) is currently used to treat motor symptoms associated with advanced Parkinson’s disease (PD) and dystonia. In particular, stimulation of the the subthalamic nucleus (STN) appears to inhibit STN hyper activity, influencing motor circuits and ultimately leading to reduce rigidity, resting tremor and hypokinesia (Benabid et al., 2009; Krack et al., 2003; Limousin et al., 1995). However, there is no clear understanding of the mechanisms of action of DBS. The STN neurons are known to display burst hyperactivity in PD patients, and one hypothesis is that DBS reduces or regulates such hyperactivity, directly or indirectly by acting on cortical projections (Gradinaru et al., 2009; Hammond et al., 2008; Montgomery and Gale, 2008). One approach to better appraise these hypotheses is to record the changes occurring in neuronal activities during DBS (Benazzouz et al., 2000; Hashimoto et al., 2003; Meissner et al., 2005). However, the neuronal signals

∗ Corresponding author at: Dept. Informatique, ESIEE-Paris, Université Paris Est, Cité Descartes-BP 99, 93162 Noisy-Le-Grand, France. Tel.: +33 145926598; fax: +33 145926699. E-mail address: [email protected] (T. Al-ani). 0165-0270/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jneumeth.2011.03.022

that can be recorded in brain structures during DBS are always contaminated by high-amplitude stimulus artefact (SA), typically spike-shaped followed by an exponential decay. We encountered the same problem by studying the impact of cortical stimulation on basal ganglia activities in a primate model of parkinsonism (Drouot et al., 2004). Generally, SA amplitude and time constant depend on several factors (O’Keeffe et al., 2001; Harding, 1991), such as experimental conditions, stimulator and recording electrode specifications, or filtering characteristics. In addition, the overlap between SA and the evoked physiological response strongly depends on the distance between the recording site and stimulus location (Harding, 1991). This distance is very short in DBS applications (as short as a few millimetres at most, Welter et al., 2004) and it affects the conduction latency period. A long, very slow exponential decay may offset SA waveform during data recording. One way to eliminate SA is to avoid amplifier saturation (Erfanian et al., 1998; Harding, 1991). There are many orders of magnitude difference between the amplitude of neuronal action potentials and the SA (SA amplitude is usually ten times more greater than action potential amplitude). The large amplitude of the SA associated with stimulation applications such as DBS, may either saturate the recording preamplifier, or contaminate the recorded neuronal activity. Therefore the dura-

136

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

tion of the usable signal between successive stimulations becomes smaller. Thus, it is important to reject or heavily attenuate the SA in neuronal recording. In this work, we propose a simple efficient method to remove SA. This method is based on the combination of two algorithms: a simple template matching algorithm of the SA and a revised version of the empirical mode decomposition algorithm (EMD) (Huang et al., 1998a) named ensemble empirical mode decomposition (EEMD) (Wu and Huang, 2009). This approach overcomes some of the limitations associated with the techniques previously used to remove SA. The main advantage of this approach is that instead of subtraction of the SA template from the original signal, the neuronal signal without stimulation artefact is reconstructed directly by eliminating the intrinsic mode functions corresponding to the SA dynamics that are generated by the EEMD algorithm. The reconstruction of the SA is made by a linear combination of a proper weighting of these modes. Then, the SA-free neuronal signal is obtained by a linear combination of a proper weighting of other modes. In this way, the dynamics of the neuronal signal are conserved. This approach should also be suitable for investigating the effects of HF-DBS (stimulation and recording at basal ganglia level) technique in the treatment of PD. 2. Background Three methods have been so far employed in the literature to eliminate the SA. 2.1. SA blanking or sample-and-hold blanking Hardware-based SA blanking methods (Black et al., 1983; Freeman, 1971; Minzly et al., 1993; Roby, 1975) are often used in conjunction with signal filtering. In this technique, a hardware circuit switches to hold mode immediately prior to stimulation to prevent the recording amplifier from detecting the SA and loosing all of the signal information during that time. At the end of the “hold” period, the circuit is switched back to sample mode and the neuronal signal again passed through to the recording amplifier. These techniques are useful in real-time applications and in preventing amplifier saturation. They are particularly useful when the SA and signal are clearly separated in time, which is not the case in HF-DBS applications. The sample-and-hold technique has largely been limited to studies using low sampling frequency electrical stimulation and those investigating long latency responses. A software versions of the sample-and-hold circuit (O’Keeffe et al., 2001; Handa et al., 1999; Keller et al., 1998) have been successfully developed for a number of experimental conditions. The use of hardware or software blanking techniques depends on the assumption that the SA and the recorded neuronal signal during stimulation are not overlapping in time and spectral domains. 2.2. Hardware/software and hybrid-based SA filtering Hardware filtering such as gain switching, slew rate limiting, or constant current/voltage switching techniques (Epstein, 1995; Parsa et al., 1998; Pozo and Delgado, 1978) and software filtering (Gnadt et al., 2003; Litvak et al., 2003; Zhang et al., 2007) algorithms were also developed to reduce the SA. All applied filters influence the quality of the recorded neuronal signal during stimulation because of the overlapping between this signal and the SA in the time and frequency domains. 2.3. SA subtraction The subtraction methods have been used to remove the SA from the recorded neuronal signal. These methods mainly include SA template subtraction (Blogg and Reid, 1990; Hashimoto et al., 2002;

Litvak et al., 2001; McGill et al., 1982; Miller et al., 1999; Wichmann, 2000), hybrid methods combining SA blanking and SA template subtraction (Montgomery et al., 2005) and removal of SA identified by local curve fitting algorithms (Wagenaar and Potter, 2002). Generating an accurate reference waveform faithfully representing the SA, is the main focus of research in subtraction techniques and is accomplished by various methods (Wagenaar and Potter, 2002; Wichmann, 2000) successfully optimising the fit of SA template to each SA event to overcome SA waveform variations. These techniques improve the effectiveness of the SA subtraction. However, Heffer and Fallon (2008) advocated that these techniques used to optimize the fit of an artefact template to each artefact event are not suitable for use within experiments investigating responses to high-rate electrical stimulation. The reason is that the computational efficiency is essential for the rapid processing of large numbers of stimulus artefact events such as HF-DBS experiments. In order to overcome the incomplete removal of the SA waveform from the recorded signal, the SA template subtraction methods assume that the individual SA waveforms match the corresponding template without a significant change in amplitude; however, this assumption cannot be made for several HFS applications (Heffer and Fallon, 2008). Miller et al. (2006) showed that the SA waveform can also vary throughout the presentation of a train of stimulus pulses. To deal with the changes in the SA waveform across stimulus intensity levels (Litvak et al., 2001; Miller et al., 2006) proposed scaling the amplitude of the SA to match the stimulus intensity. However, because the relationship between stimulus intensity and SA waveform is complex and often unknown, this technique can result in residual SA remaining in the recorded signal (Heffer and Fallon, 2008; Litvak et al., 2001). Heffer and Fallon (2008) proposed an effective sample-and-interpolate SA rejection technique suitable for high sampling rate where the recorded signal is digitized to capture the shape of the neuronal response and the timing of the SA events. This technique successfully considers that the SA and the recorded neuronal signal during stimulation are overlapping in time and spectral domains. However, as indicated by Heffer and Fallon (2008), this technique requires that the artefact events be identifiable and that the artefact duration remains less than both the inter-stimulus interval and the time course of the neuronal response. The EMD (Huang et al., 1998a) and its revised version, the EEMD (Wu and Huang, 2009) are proposed as relatively novel techniques which overcome the limitations of the Wigner–Ville distribution (Ville, 1948), and the short-time Fourier transform for processing nonlinear and nonstationary signals. Flandrin et al. (2004a) reported that EMD behaves as a “wavelet-like” filter bank. The advantage of EMD, with respect to wavelet transform (Mallat, 1989), is that it is a data driven approach, i.e. one does not need to define a mother wavelet beforehand. The EMD is used successfully in biomedical signal analysis (Balocchi et al., 2004; Blanco-Velasco et al., 2008; Echevarría et al., 2001; Flandrin et al., 2004b; Huang et al., 1998b; Liang et al., 2000, 2005; Salisbury and Sun, 2004). The EMD is applied with success to signal noise denoising and to reduce SA (Blanco-Velasco et al., 2008; Lieu et al., 2008).

3. Methods and results 3.1. Record procedure The SA removal algorithm was developed to help in the analysis of electrical neuronal signals, which were recorded during surgical implantation of DBS electrodes in the STNs in a PD patient who was previously treated by implanted motor cortex stimulation (MCS). A stainless steel microelectrode (Medtronic, Minneapolis, USA) was inserted through a burr hole into a guide cannula according to

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

137

Fig. 1. Diagram of the experimental procedure. A microelectrode was implanted in the right STN and another one in the left STN. Neuronal STN activity was recorded, while the right motor cortex was stimulated using a previously implanted epidural electrode.

a stereotactic trajectory of which end point was chosen to be the infero-lateral part of the STN. The microelectrode was advanced using a manual microdrive device and spontaneous neuronal activity was recorded along the trajectory with a Leadpoint analysis station (Medtronic). The signal was amplified, bandpass filtered at 300–5000 Hz and sampled at 24 kHz. Then, the electrode was placed in the zone of maximal STN hyperactivity, according to the typical electrophysiological features of STN cell activity (Benazzouz et al., 2002; Hutchison et al., 1998; Pralong et al., 2002). At this optimal location, STN activity was recorded for bins of 30 s, before, during, and after that the previously implanted epidural MCS was switched ‘on’, for a total duration of 5 min. MCS parameters were as follows: frequency 40 Hz, pulse duration 150 ␮s, voltage 3 V. The procedure was performed twice, one for recordings in the right STN and the other for recordings in the left STN (Fig. 1). In this work, we were interested only in the 30-s records obtained during the stimulation. Five records from right STN and four records from left STN were used in this study. Each 30-s record is composed of a number of segments as seen in Fig. 2. Each segment (about 12 ms length) is composed of two consecutive time series signals: an SA-free neuronal signal (about 10 ms length) followed by one SA (about 2 ms length). Whatever the parameters of stimulation, all the artefacts were always within the working range of the amplifier used in our recording system. 3.2. EMD and EEMD approaches In this section, we give a brief introduction to the basic principles of EMD and the EEMD techniques used in our SA removal algorithm. 3.2.1. EMD approach The traditional EMD was recently proposed (Huang et al., 1998a) as an adaptive time-frequency data analysis method. It is defined by an algorithm based on an empirical framework. In most cases, the studies (performance, analysis, etc.) carried out on the EMD are done with extensive digital simulations in controlled conditions

(Huang et al., 1998a). Despite the lack of theoretical formalism, this algorithm showed its capacity to analyse the signals. Using a new formulation for EMD based on constrained optimization, the results of Meignen and Perrier (2007) were very similar to those obtained with the traditional EMD algorithm. The basic EMD is defined by a process called sifting to break down any multimodal signal to a sum of basis components called intrinsic mode functions (IMFs). The IMFs are zero-mean AM-FM signals which must satisfy two conditions: the first one is that the number of extrema and that of zero-crossing must differ at most by one; the second one is that the mean value between the upper and lower envelopes are equal to zero at any point. Conceptually, the establishment of this method is quite simple: one needs to consider a signal at its local oscillation level, remove the fastest oscillation and iterate the process on the residue considered as a new signal. At the end of the sifting processes, a given signal x(t) can be written as a sum of a finite number of IMFs, Im (t), m = 1, 2, . . ., M, and a final residue rM (t): x(t) =

M 

Im (t) + rM (t).

(1)

m=1

The decomposition is stopped at step M, if either the residue rM (t) is a mono-component signal or has only 2 extrema (Huang et al., 1998a). The stopping criterion must be set to ensure that the obtained signal satisfies the properties of an IMF while limiting the number of iterations. For more details about the different steps of the sifting process for the calculation of the IMFi as well as the stopping criterion definition (see Huang et al., 1998a). Since the decomposition into IMFs is based on the local characteristic time scale of the data, it applies to nonlinear and non-stationary processes. 3.2.2. EEMD approach The EEMD is a revised version of the EMD that gives a substantial improvement over the original EMD. One of the major drawbacks

138

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

Fig. 2. Template and EMD/EEMD based SA removing algorithm.

of EMD is mode mixing; it is a consequence of signal intermittency. Mode mixing is defined as a single IMF either consisting of signals of widely disparate scales, or a signal of a similar scale residing in different IMF components.

Mode mixing could not only cause serious aliasing in the time-frequency distribution, but also make the physical meaning of individual IMF unclear. To alleviate this drawback, a new noise-assisted data analysis (NADA) method,

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

the Ensemble EMD (EEMD) is proposed by Wu and Huang (2009). The EEMD defines the true IMF as the average over an ensemble of tests; each test is the EMD of the original signal added to white noise. The added white noise is treated as the possible random noise that would be encountered in the measurement process. Under such conditions, the kth “artificial” observation will be xk (t) = x(t) + wk (t),

(2)

This procedure produces a collection of white noises which cancel each other; therefore, only the real components can survive and persist in the final average. The amplitude of white noises must force the ensemble to get all possible solutions: the noise makes the various components reside in the corresponding IMF dictated by the EMD filter banks and the significant physical sense of IMF. The EEMD is a time consuming procedure due to the high number of tests (typically 100). 3.3. Template and EEMD-based SA removal algorithm Our SA removal algorithm was developed using a simple direct template extraction of the artefact and the EEMD program developed by Wu and Huang (2009). The details of this algorithm are shown in Fig. 2. The algorithm may be summarised by the following steps: 1. Raw input signal (Fig. 2(a)): The raw input signal is recorded during stimulation. It is composed of successive segments of neuronal activity and a high amplitude SA. In order to keep the frequencies of interest, this signal must be analysed by a common spectral analysis method. 2. Band-pass filtering (Fig. 2(b)): Based on the spectral analysis, the raw input signal is first pre-filtered using a band-pass filter to keep frequencies of interest between a minimum frequency Fmin and a maximum frequency Fmax . In our work, these frequencies were 300 Hz and 5000 Hz respectively, which correspond to the usual frequency bandpass for recording STN cell activity (e.g., Hashimoto et al., 2002; Pralong et al., 2002). A stream of input signals issued from one session is read on-line. 3. On-line SA template and pre-filtered input segment extraction (Fig. 2(c)): As discussed in Section 1, the morphology of the SA depends on several factors and especially the nature of the stimulation experiment. In this work, we are interested on SA produced during our MCS neurosurgical technique. The produced SA has the same morphology everywhere in the analysed stimulation signals. Thus, the simple template extraction method presented in this paragraph is closely related to this type of SA. The following steps summarise the current on-line SA detection and extraction. (a) start by a sample index i0 , explore the pre-filtered signal by incrementing the sample index until a value xp (ip ) greater than or equal to a given positive threshold thp is found. (b) start from the sample index ip , decrement the sample index until a left-side minimum value xminL (iminL ) less than or equal to a given left-side negative threshold (thnL ) is found. (c) start from the sample index iminL , decrement the sample index until a left-side zero-crossing is found. The associated index izcL is considered to be the starting index of the current SA. (d) start from the sample index ip , increment the sample index, stock the index izcR = izcR1 corresponding to the first rightside zero-crossing. (e) continue incrementing the sample index until a right-side minimum value xminR (iminR ) is found. • if this value is less than or equal to a given right-side negative threshold (thnR ) then start from the sample index

139

(iminR ), increment the sample index until a right-side zerocrossing is found. The associated index (izcR = izcR2 ) is considered to be the ending index of the current SA. • else, take the index izcR1 as the ending index of the current SA. Thus, the last extracted segment is given by the part of the prefiltered input signal defined by the sample indexes i0 and izcR and the last extracted template is defined by the sample indexes izcL and izcR of this segment. The SA template in the case of DBS may be extracted in the same manner with a little modification of the above method. Considering separately the fluctuations of periodograms for the two parts of a given segment, i.e., neuronal activity before the SA and during the SA, the mean power spectral density may be then calculated. Fig. 3 shows an example of the Welch power spectral densities of these two parts in the signal given in Fig. 2(a). In this example, we noted that there is a significant overlap between the spectral contents of the neuronal signal and the SA. This overlap makes the conventional filters ineffective for eliminating the SA. 4. Cross-correlation (Fig. 2(d)): Check the consistency of the extracted template by calculating the cross-correlation between the last extracted template and the current segment. 5. Window function (Fig. 2(e)): The cross-correlation must be greater than a given threshold thc at the SA zone between izcL and izcR . Once this correlation is satisfied, a window function is constructed to further preserve the IMFs in the SA zone and it is calculated as follows:

 W (t) =

1 if t(izcL ) ≤ t ≤ t(izcR ) and C > thc , 0, otherwise.

6. EMD/EEMD decomposition (Fig. 2(f)): Apply the EEMD algorithm on the extracted input segment. This segment is then decomposed into several IMFs, as mentioned in Sections 3.2.1 and 3.2.2, with increasing indexes from high to low frequency where the highest index correspond to residual (over all trend). Using the EEMD approach, one can separate scales naturally without any a priori subjective criterion selection (Wu and Huang, 2009). In this work, we proposed an on-line artefact processing that may be applied on short segments of the neuronal signal during stimulation. A segment is composed of neuronal activity followed by an SA (Fig. 2(c and f)). As seen in Fig. 2(f), the frequency content of each individual IMF decreases as the order of IMF increases. Note that the high order IMFs are mainly due to the SA, which has strong high-frequency components. This observation can be used to delineate the SA. To delineate the IMFs of the SA and of the SA-free neuronal activity, frequency spectrum analysis by periodogram may be also adopted. Considering the fluctuations of periodogram for the given pre-filtered signal, the mean power spectral density may be calculated. Fig. 4 shows an example of the Welch power spectral densities of the IMFs given in Fig. 2(f). 7. SA and SA-free neuronal signal reconstitution (Fig. 2(g–i)): The  may be obtained as the weighted sum of reconstructed SA, SA, the IMFs on the SA zone defined by the window function:

 = SA

M 

a(m)W (t)Im (t),

(3)

m=1

where M is the number of IMFs, Im (t) is the IMF corresponding to the index m, and a(m), 0 ≤ a(m) ≤ 1, is a weight coefficient chosen to decay gradually. In our work, the values of the weight coefficient

140

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

Welch power spectral density (dB/Hz) of the analysed pre−SA neuronal activity 0.1 0.05 0

0

2000

4000

6000

8000

10000

12000

10000

12000

Frequency (Hz) Welch power spectral density (dB/Hz) of the SA signal

10 5 0

0

2000

4000

6000

8000

Frequency (Hz) Amplitude (mV)

The complete analysed frame: pre−stimulus neural signal and the stimulation artefact 500 0 −500

0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

Time (s) Fig. 3. Welch power spectral densities (dB/Hz) of the neuronal activity (Top) and of the SA (middle) extracted from the segment shown in the bottom of the figure.

(a) Normalised root mean squared error (NRMSE) The RMSE also evaluates the accuracy which is defined by

are chosen as follows: a(m)

10 − m , if m = 1 10 K −m = otherwise. K =

  N  (y(i) − yˆ (i))2 RMSE =  .

(4)

This choice is inspired by the fact that the SA zone spreads with increasing order of IMF. The SA-free neuronal signal, xˆ r , is reconstructed by xˆ r =

M 

 Im (t) − SA.

(5)

m=1

ˆ using the current extracted The coefficient K is estimated (K) segment by optimising one of the following standard quantitative mathematical errors formulations measures: Let y, the true sub-segment (SA-free or SA) of the input segment and yˆ the reconstructed sub-segment (SA-free or SA) using our algorithm.

(6)

N

i=1

The normalised RMSE (NRMSE) is given by NRMSE = RMSE/(ymax − ymin ),

(7)

(b) Percent root mean square difference (PRD) The PRD, given in percentage, indicates reconstruction fidelity by point-wise comparisons with the original data which is defined by



N

PRD = 100 ×

(y(i) − yˆ (i))

i=1 N

i=1

¯ 2 (y − y)

2

,

(8)

where ymax , ymin and y¯ are the maximum, the minimum and the mean of y respectively.

Welch power spectral densities (dB/Hz) of IMF1 to IMF7, of residue (red/dashed line), and of pre−filtered signal (blue/full line) 0.04 0.02 0 0 0.1 0.05 0 0 0.1 0.05 0 0 0.1 0.05 0 0 0.04 0.02 0 0 0.04 0.02 0 0 0.1 0.05 0 0 0.2 0.1 0 0 0.4 0.2 0 0

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

1000

2000

3000

4000

5000

6000

Frequency (Hz) Fig. 4. From top to down: Welch power spectral densities (dB/Hz) of IMF1 to IMF7 (red/dotted line), of residue (red/dotted line), and of the pre-filtered input segment (blue/full line) given in Fig. 2(f). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146 Signal with artefact (blue/full line) reconstructed artefact by emd (red/dashed line)

200

200

100

100

Amplitude (mV)

Amplitude (mV)

Signal with artefact (blue/full line) and the reconstructed artefact by emd (red/dashed line)

141

0

−100

−200

0

−100

−200

−300 0

0.002

0.004

0.006

0.008

0.01

0.012

−300 0

0.014

0.002

0.004

0.006

Time (s)

1.4

1.2

1.2

1 0.8 0.6 0.4 0.2 0

2000

4000

6000

0.01

0.012

0.014

Power densities: pre−filtered (blue/full line) and clean reconstructed signal (red/dashed line)

1.4

Power/Frequency (dB/Hz)

Power/Frequency (dB/Hz)

Power densities: SA template (blue/full line) and reconstructed artefact (red/dashed line)

0

0.008

Time (s)

8000

10000

1 0.8 0.6 0.4 0.2 0

12000

0

2000

4000

6000

8000

10000

12000

Frequency (Hz)

Frequency (Hz) Signal with artefact (blue/full line) reconstructed artefact by emd (red/dashed line)

Amplitude (mV)

40 20 0 −20 −40 −60 0.0095

0.01

0.0105

0.011

0.0115

0.012

0.0125

Time (s) Fig. 5. Top-left: the true signal with artefact (blue/full line) and the reconstructed SA by emd, SA (red/dashed line). Top-right: the true signal with artefact (blue/full line) and

the reconstructed SA-free neuronal signal, xˆ (red/dashed line). Middle-left: Welch power spectral densities of the true signal with SA (blue/full line) and of SA (red/dashed line). Middle-right: Welch power spectral densities of the true signal with SA (blue/full line) and of xˆ (red/dashed line). Down: Zoom on xˆ . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

(c) Spectral relative error (SRE) The SRE, given in percentage, quantifies the performance of the EMD/EEMD to reconstruct the artefact and the neuronal signals.

Fmax SRE =

Ryˆy =

N

N

(Py f ) − Pˆ yˆ (f ))2

f =Fmin Fmax P 2 (f ) f =Fmin y



The Ryˆy (0 ≥ Ryˆy ≤ 1), assesses the similarity between the realistic and reconstructed signals.

,

(9)

N i=1

i=1

y(i)ˆy(i) −

N

y2 (i) − (

i=1

y(i))

N N y(i) yˆ (i) i=1

i=1 N 2 N 2 N

i=1

yˆ (i) − (

i=1

.

(10)

2

yˆ (i))

(e) Signal to error ratio (SER) The SER gives the error measures which is simply defined by

N

SER = 10 log where Py (f) is the spectral density of the given signal segment y, Pˆ yˆ (f ) is the spectral density of the reconstructed signal yˆ and f ranges from Fmin to Fmax . (d) Pearson product-moment correlation coefficient (Ryˆy )

N

N i=1

i=1

y2 (i)

(y(i) − yˆ (i))

2

.

(11)

(f) Peak signal-to-noise ratio (PSNR) The PSNR used to describe objectively the quality of the reconstructed data. PSNR is the ratio between the maximum possible

142

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

Table 1 Statistics of the the six performance measures (PM) using EMD and EEMD. PM

Statistics of the PM

EMD

NRMS

Mean for the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

0.1283 0.0853 0.2145 0.0739

0.0554 0.0258 0.0781 0.0071

PRD

Mean of the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

0.9243% 0.4487% 0.7694% 0.2104%

0.4233% 0.1875% 0.3461% 0.0182%

SRE

Mean for the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

0.3676 0.6164 1.3749 4.0998

0.0035 0.0030 0.0533 0.0191

R

Mean for the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

0.5698 0.2239 0.7024 0.1660

0.8884 0.0892 0.9607 0.0046

SER

Mean for the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

1.3520 dB 3.0674 dB 2.1856 dB 4.0433 dB

8.6478 dB 4.7388 dB 10.3660 dB 0.3281 dB

PSNR

Mean for the post-SA neuronal signal Std for the post-SA neuronal signal Mean for the SA Std for the SA

15.5604 dB 2.0805 dB 11.2941 dB 3.4984 dB

21.6318 dB 4.0208 dB 19.1106 dB 0.7497 dB

power of the given signal segment and the power of corrupting noise that affects the fidelity of its representation. A higher PSNR value indicates that the reconstruction is of high quality. Acceptable values are considered to be about 20–25 dB.

N PSNR = 10 log

Ni=1 i=1

¯ (y(i) − y)

2

(y(i) − yˆ (i))

2

.

(12)

A small RMSE, PRD, SRE and with a large ryˆy , SER and PSNR indicate good performance in SA and SA-free neuronal signals. The use of these measures imply that the processed segment (clean neuronal signal and SA) has no baseline fluctuations as the variance of this segment stays almost constant. Algorithm 1 summarises the selection of this parameter for the current pre-filtered input segment: Algorithm 1. Estimation of the coefficient K for the given segment. Select one performance measure (PM) among the performance measures given above. Input data: window function (W(t)); extracted pre-filtered segment (segm); IMFs matrix (M × length(segm)); extracted template’s length (Ltempl ); maximum iteration number (Kmax ) (e.g., 1000). Initialisation: m = 1 (first IMF), K1 = 10 or a(1) = 0.9 (using Eq. (4)) for K = 100 to Kmax do Calculate a(m) for m = 2 to m = M (using Eq. (4)) Estimate xˆ r and SA (using Eq. (3) and Eq. (5)) Calculate and stock the selected performance measure PM(K). end for Select Kˆ such that Kˆ = arg min (or max)(PM(K)) Use Kˆ to estimate the optimal values of xˆ r and SA (using Eq. (3), Eq. (4) and Eq. (5))

As an example, Fig. 5 shows respectively a reconstructed SA,  and a reconstructed SA-free neuronal signal, xˆ with their associSA, ated power spectral densities using the above approach. This figure shows the reliability of our approach.

EEMD

 (pre- or post-SA).  In the reconstructed neuronal signal, xˆ outside SA same manner, we compared the true extracted SA template with  The minimal (or maximal) values of each the reconstructed SA (SA). performance measure calculated by Algorithm 1 were stocked for each segment of the current processed post-filtered signal. Table 1 shows a comparative results between emd and eemd decompositions using these performance measures applied on the same input signal. The input signal of 1 s length used to obtain these results is composed of 81 segments (pre-SA neuronal signal followed by one SA signal). The results show that all the selected performance measures give good statistical results for the two types of decompositions. However, we recommend the use of eemd decomposition because it gives better performance measures in the presence of noise and probably hidden mode mixing. All the performance measures gave almost the same results, particularly NRMS performance measure. Fig. 6 gives the results in this case. A Box and whisker plot (McGill et al., 1978) is presented to show the distribution. The expanded traces shown in Fig. 6 demonstrate that our method of artefact removal does not degrade the signalto-noise ratio of the data and does not alter parts of the signal in which artefacts were not found. The boxes have lines at the lower quartile, median, and upper quartile values. The whiskers are lines extending from each end of the box to the extreme data values not defined as outliers. Points lying beyond 1.5 times the interquartile range from the median are shown as outliers (+). 4. Discussion

3.4. Performances and validation

Stimulus artefacts (SA) often contaminate recorded neuronal signals and make the analysis of the neuronal response in DBS applications difficult. The technique used in our approach removes, on-line, the SA from contaminated neuronal signal by requiring the data to pass some number of application-related threshold points. We may consider the following remarks:

To evaluate the effectiveness of our approach, we compared the true neuronal signal, x, outside the SA (pre- or post-SA) with the

1. We used three threshold points. A positive threshold thp , on the positive peak, corresponds to about the half of this peak. Two

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

NRMSE of the reconstructed pre−SA neuronal signal

143

NRMSE of the reconstructed pre−SA neuronal signal 0.2 0.15

0.4

NRMSE

NRMSE

0.5

0.3 0.2

0.1 0.05

0.1 1

1

NRMSE of the reconstructed SA signal 0.1

0.3

NRMSE

NRMSE

0.4

NRMSE of the reconstructed SA signal

0.2

0.08

0.1 0.06

1

1

Pr−filtered signal with SA (blue/full line) and post−filtered SA−free signal (red/dashed line)

300

300

200

200

Amplitude (mV)

Amplitude (mV)

Pr−filtered signal with SA (blue/full line) and post−filtered SA−free signal (red/dashed line)

100

0

−100

0 −100 −200

−200 −300

100

0

0.01

0.02

0.03

0.04

0.05

0.06

0.07

0.08

0.09

−300

0.1

0

Time (s)

100

100

Amplitude (mV)

Amplitude (mV)

200

0.04

0.05

0.06

0.07

0.08

0.09

0.1

0

−100

−200

−200

−300

0.03

Pr−filtered signal with SA (blue/full line) and post−filtered SA−free signal (red/dashed line)

200

−100

0.02

Time (s)

Pr−filtered signal with SA (blue/full line) and post−filtered SA−free signal (red/dashed line)

0

0.01

0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Time (s)

−300

0

0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01

Time (s)

Fig. 6. NRMS performance measure. Left column: results using EMD approach. Right column: results using EEMD approach. For each column and from top to down: Box and whisker plot, a part (100 ms) of the pr-filtered signal with SA (blue/full line) and post-filtered SA-free signal (red/dashed line), zoom on the first segment of this signal and an expanded traces of the SA-free pr-filtered signal and the reconstructed SA-free post-filtered signal. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

negative thresholds: thnL and thnR , on the left-side and rightside negative peaks respectively. thnL was set to about the half of the left-side negative peak, and thnR was set to about 10% of the left-side negative peak. The number and the values of these thresholds may be set manually and easily in any MCS or DBS applications where the SA morphology is stable during all stimulation periods. To make our on-line SA extraction method more flexible to work with other type of recorded stim-

ulus contaminated electroneurophysiologic data, the reader may combine the approach given by O’Keeffe et al. (2001) with our approach. 2. In order to maximise the efficiency of the algorithm to remove only the SA, we applied dynamically a threshold-related error control. The program calculates the mean of the upper and lower envelops of the pre-SA clean neuronal wave and compare, at the moment where the first threshold is detected, the three thresh-

144

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

300

Amplitude (mV)

200 100 0 −100 −200 −300 −400

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

1.8

2

Time (s) 300

Amplitude (mV)

200 100 0 −100 −200 −300 −400

0

0.2

0.4

0.6

0.8

1

Time (s) Fig. 7. A part (2 s: 162 segments (SA-free signal and SA)) of about 21,600 segments used to test the robustness of the reconstructed SA-free signal. Up: The true signal with SA (blue/full line) and the reconstructed SA (red/dashed line). Down: The true signal with SA (blue/full line) and the reconstructed SA-free signal (red/dashed line). A zoom on one reconstructed SA and one reconstructed SA-free signal are shown in Fig. 5. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

old values to the means of these envelops. These means must never exceed thp , thnL and thnR . 3. The performance of the artefact removal technique is not affected by the firing rate. HFS does not change the shape of the artefact, simply the number of SA will be increased. The procedure for extraction of template remains unchanged. In addition, the sifting process of the EMD algorithm is independent of the frequency range of the input signal. Finally, the method does not alter parts of the signal in which artefacts were not found, as shown in Fig. 6. 4. Our program was developed using MATLAB (Mathworks). However, to use it in on-line, it may be easily implemented in C language and compiled as a standalone executable file. A realtime implementation of our algorithm may be achieved using a computer with multi-core processors or a real-time digital signal processing (DSP). 5. The EMD is a better option if processing time is a limitation (using moderate computing system) and no mode mixing phenomena (Wu and Huang, 2009) is present in the input signal (which is the case in our analysed segments). In this case, the noise must be filtered before any decomposition. The EMD approach may be efficiently used for signal de-noising (e.g., Gang et al., 2006; Kopsinis and McLaughlin, 2009). In our work, the input signals contained a very low amount of noise. However, the de-noising step was realised at the same time as the SA processing by removing the high frequency components given by the higher order IMFs with minimum signal distortion.

When using EEMD, the final EMD is estimated by averaging numerous EMD runs with the addition of noise. This technique was an advancement introduced by Wu and Huang (2009) to try increasing the robustness of EMD in the presence of noise and alleviate some of the common problems of EMD such as mode mixing. Niazy et al. (2009) show that EEMD, in addition to slightly increasing the accuracy of the EMD output, substantially increases the robustness of the results and the confidence in the decomposition. In our work, the very low amount of noise present in the input signals was removed directly using the EEMD. The EEMD is more time consuming than EMD when using Matlab (Mathworks). However, the implementation of our algorithm in C language or on DSP for on-line processing may alleviate this problem. 6. In our experiment, the amplitude and shape of the SA were substantially the same during all stimulation periods. As described in Section 3.3, our on-line SA template detection algorithm extracts directly the SA in the current processed segment (SA free neuronal signal and SA) and considers it as a template without using any general (or model-based) template matching. Hence, the application of our algorithm to all recorded data gave a false error rate of 0% for the SA-free neuronal signal as well as for the SA. Thus, in order to evaluate the robustness of our approach, the respective number of SA and SA-free neuronal signal was compared between the true and the reconstructed signals. The comparisons were done on the basis of nine signals (nine 30-s records) recorded in a single patent with a total number of about

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

21,600 segments. The accuracy of the obtained results demonstrates the robustness of our proposed algorithm (see Fig. 7 as an example). However, for experiments in which the SA could change in amplitude and/or shape; a fixed or dynamic template should be needed. Hence, only the SA template extraction procedure described in step 3 of Section 3.3 should be modified. In this case, the robustness evaluation of our algorithm should be evaluated by estimating the expected false positive and false negative error rates. 7. The main advantages of our algorithm may be summarized as follows: (a) Our algorithm avoids the use of analogue or digital filtering techniques that need a priori information on signal’s frequency. (b) As discussed previously, for DBS and MCS experiments, where the SA does not change significantly in amplitude and in shape, the SA template is extracted directly as the true on-line processed SA. This technique makes the SA removal process highly robust. (c) Our technique does not degrade the signal-to-noise ratio of the data and does not alter parts of the signal in which artefacts were not found (as discussed previously). (d) The reconstruction procedure of the SA and SA-free signals does not need any a priori information on signal’s frequency and stimulus timing. The selection of the IMFs corresponding to these signals is based on an optimisation step in order to minimise some error criterion between the true and the reconstructed signals (as described in Section 3.3). The main drawback of our algorithm may be summarised as follows: (a) The SA template extraction step is application-dependent. However, in order to include other stimulation experiments with different SA waveforms this step may be adapted independently of the other steps of our algorithm. (b) The reconstruction procedure of the SA and SA-free signals is based on the EMD or the EEMD and may be time consuming. However, in order to accelerate the computational process for on-line or real (or near-real) time processing, the algorithm may be implemented in C language or on a DSP (as discussed previously). This could offer more flexibility throughout the stimulation experiment. 5. Conclusions In this work, we demonstrated efficiency of the direct template extraction and the EMD/EEMD for SA removal as applied to neuronal activity recorded in the STN during MCS in a PD patient. Our approach may also be easily adapted and applied to the case of any basal ganglia recording while HF-DBS is switched ‘on’. In the case of HF-DBS, stimulus artefact should be removed for frequency of 100 Hz or higher (typically 130 Hz) and not only for 40 Hz as in this study. However, higher frequencies of stimulation do not alter the performance of our method as explained above. Application of this approach enables one to estimate the firing rate and further response analysis during stimulation periods, in order to investigate more clearly the physiological mechanisms of action of any implanted neurostimulation technique into the brain. Acknowledgments The authors gratefully thank Drs Bechir Jarraya, Hélène Lepetit, Jean-Marc Gurruchaga, Gilles Fénelon, Xavier Drouot, and Philippe Rémy.

145

References Balocchi P, Menicucci D, Santarcangelo E, Sebastiani L, Gemignani A, Ghelarducci B, et al. Deriving the respiratory sinus arrhythmia from the heartbeat time series using empirical mode decomposition. Chaos Solitons Fractals 2004;20:171–7. Benabid AL, Chabardes S, Mitrofanis J, Pollak P. Deep brain stimulation of the subthalamic nucleus for the treatment of Parkinson’s disease. Lancet Neurol 2009;8:67–81. Benazzouz A, Gao DM, Ni ZG, Piallat B, Bouali-Benazzouz R, Benabid AL. Effect of highfrequency stimulation of the subthalamic nucleus on the neuronal activities of the substantia nigra pars reticulata and ventrolateral nucleus of the thalamus in the rat. Neuroscience 2000:289–95. Benazzouz A, Breit S, Koudsie A, Pollak P, Krack P, Benabid AL. Intraoperative microrecordings of the subthalamic nucleus in Parkinson’s disease. Mov Disord 2002;17(Suppl. 3):S145–9. Black RC, Clark GM, ÓLeary SJ, Walters C. Intracochlear electrical stimulation of normal and deaf cats investigated using brainstem response audiometry. Acta Otolaryngol Suppl 1983;399:5–17. Blanco-Velasco M, Weng B, Barner KE. ECG signal denoising and baseline wander correction based on the empirical mode decomposition. Comput Biol Med 2008;38:1–13. Blogg T, Reid WD. A digital technique for stimulus artifact reduction. Electroencephalogr Clin Neurophysiol 1990;76:557–61. Drouot X, Oshino S, Jarraya B, Besret L, Kishima H, Remy P, et al. Functional recovery in a primate model of Parkinson’s disease following motor cortex stimulation. Neuron 2004;44:769–78. Epstein CM. A simple artifact-rejection preamplifier for clinical neurophysiology. Am EEG Technol 1995;35:64–71. Echevarría JC, Crowe JA, Woolfso MS, Hayes-Gill BR. Application of empirical mode decomposition to heart rate variability analysis. Med Biol Eng Comput 2001;39:471–9. Erfanian A, Chizeck HJ, Hashemi RM. Using evoked EMG as a synthetic force sensor of isometric electrically stimulated muscle. IEEE Trans Biomed Eng 1998;45:188–202. Flandrin P, Rilling G, Goncalves P. Empirical mode decomposition as a filter bank. IEEE Signal Process Lett 2004a;11:112–4. Flandrin P, Goncalves P, Rilling G. Detrending and denoising with empirical mode decomposition. In: Proceedings of XIIth EUSIPCO, 2004; 2004b. Freeman JA. An electronic stimulus artifact suppressor. Electroencephalogr Clin Neurophysiol 1971;31:170–2. Gnadt J, Echols S, Yildirim A, Zhang H, Paul K. Spectral cancellation of microstimulation artifact for simultaneous neural recording in situ. IEEE Trans Biomed Eng 2003;50:1129–35. Gradinaru V, Mogri M, Thompson KR, Henderson JM, Deisseroth K. Optical deconstruction of parkinsonian neural circuitry. Science 2009;324:354–9. Gang L, Xinwei W, Lihua S. On EMD based adaptive de-noising method for signal processing. J Adv Sci 2006;18(1/2):144–51. Hammond C, Ammari R, Bioulac B, Garcia L. Latest view on the mechanism of action of deep brain stimulation. Mov Disord 2008;23:2111–21. Handa T, Takahashi H, Saito C, Handa Y, Ichie M, Kameyama J, et al. Development of an FES system controlled by EMG signals. In: Proceedings of the IEEE 20th Conference of Engineering in Medicine and Biology Society; 1999. Harding GW. A method for eliminating the stimulus artifact from digital recordings of the direct cortical response. Comp Biomed Res 1991;24:183–95. Hashimoto T, Elder C, Vitek J. A template subtraction method for stimulus artifact removal in high-frequency deep brain stimulation. J Neurosci Methods 2002;113:181–6. Hashimoto T, Elder CM, Okun MS, Patrick SK, Vitek JL. Stimulation of the subthalamic nucleus changes the firing pattern of pallidal neurons. J Neurosci 2003;23:1916–23. Heffer LF, Fallon JB. A novel stimulus artifact removal technique for high-rate electrical stimulation. J Neurosci Methods 2008;170:277–84. Hutchison WD, Allan RJ, Opitz H, Levy R, Dostrovsky JO, Lang AE, et al. Neurophysiological identification of the subthalamic nucleus in surgery for Parkinson’s disease. Ann Neurol 1998;44:622–8. Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc R Soc London 1998a;454:903–95. Huang W, Shen Z, Huang NE, Fung YC. Engineering analysis of biological variables: an example of blood pressure over 1 day. Proc Natl Acad Sci USA 1998b;95:4816–21. Keller T, Curt A, Popovic MR, Dietz V, Signer A. Grasping in high lesioned tetraplegic subjects using the EMG controlled neuroprosthesis. J Neurol Rehabil 1998;10:251–5. Kopsinis Y, McLaughlin S. Development of EMD-based denoising methods inspired by wavelet thresholding. IEEE Trans Signal Process 2009;57(April (4)):1351–62. Krack P, Batir A, Van Blercom N, Chabardes S, Fraix V, Ardouin C, et al. Fiveyear follow-up of bilateral stimulation of the subthalamic nucleus in advanced Parkinson’s disease. N Engl J Med 2003;349:1925–34. Liang H, Lin Z, McCallum R. Artifact reduction in electrogastrogram based on empirical mode decomposition method. Med Biol Eng Comput 2000;38:35–41. Liang H, Lin QH, Chen JDZ. Application of the empirical mode decomposition to the analysis of esophageal manometric data in gastroesophageal reflux disease. IEEE Trans Biomed Eng 2005;52:1692–701. Lieu S, He Q, Gao RX, Freedson P\. Empirical mode decomposition applied to tissue artifact removal from respiratory signal. In: International Conference of the IEEE on Engineering in Medicine and Biology Society; 2008. p. 3624–7.

146

T. Al-ani et al. / Journal of Neuroscience Methods 198 (2011) 135–146

Limousin P, Pollak P, Benazzouz A, Hoffmann D, Le Bas JF, Broussolle E, et al. Effect of parkinsonian signs and symptoms of bilateral subthalamic nucleus stimulation. Lancet 1995;345:91–5. Litvak L, Delgutte B, Eddington D. Auditory nerve fiber responses to electric stimulation: modulated and unmodulated pulse trains. J Acoust Soc Am 2001;110:368–79. Litvak L, Smith Z, Delgutte B, Eddington D. Desynchronization of electrically evoked auditory-nerve activity by high-frequency pulse trains of long duration. J Acoust Soc Am 2003;114:2066–78. Mallat SG. A theory for multiresolution signal decomposition: the 619 wavelet representation. IEEE Trans Pattern Anal Mach Intell 1989;11(July):674–93. McGill R, Tukey JW, Larsen WA. Variations of boxplots. Am Stat 1978;32: 12–6. McGill KC, Cummins KL, Dorfman LJ, Berlizot BB. On the nature and elimination of stimulus artifact in nerve signals evoked and recorded using surface electrodes. IEEE Trans Biomed Eng 1982;29:129–37. Meignen S, Perrier V. A new formulation for empirical mode decomposition based on constrained optimization. IEEE Signal Process Lett 2007;14:932–5. Meissner W, Leblois A, Hansel D, Bioulac B, Gross CE, Benazzouz A, et al. Subthalamic high frequency stimulation resets subthalamic firing and reduces abnormal oscillations. Brain 2005;128:2372–82. Minzly J, Mizrahi J, Hakim N, Liberson A. Stimulus artefact suppressor for EMG recording during FES by a constant-current stimulator. Med Biol Eng Comput 1993;31:72–5. Miller C, Abbas P, Robinson B, Rubinstein J, Matsuoka A. Electrically evoked single fiber action potentials from cat: responses to monopolar, monophasic stimulation. Hear Res 1999;130:197–218. Miller C, Abbas P, Robinson B, Nourski K, Zhang F, Jeng F. Electrical excitation of the acoustically sensitive auditory nerve: single-fiber responses to electric pulse trains. J Assoc Res Otolaryngol 2006;7:195–210. Montgomery E, Gale J, Huang H. Methods for isolating extracellular action potentials and removing stimulus artifacts from microelectrode recordings of neurons requiring minimal operator intervention. J Neurosci Methods 2005;144: 107–25.

Montgomery E, Gale J. Mechanisms of action of deep brain stimulation (DBS). J Neurosci Biobehav Rev 2008;32:388–407. Niazy RK, Beckmann CF, Brady JM, Smith SM. Performance evaluation of ensemble empirical mode decomposition. Adv Adapt Data Anal 2009;1(2):231–42. O’Keeffe DT, Lyons GM, Donnelly AE, Byrne CA. Stimulus artifact removal using a software-based two-stage peak detection algorithm. J Neurosci Methods 2001;109:137–45. Parsa V, Parker P, Scott R. Convergence characteristics of two algorithms in non-linear stimulus artefact cancellation for electrically evoked potential enhancement. Med Biol Eng Comput 1998;36:202–14. Pozo FD, Delgado JMR. Hybrid stimulator for chronic experiments. IEEE Trans Biomed Eng 1978;25:92–4. Pralong E, Ghika J, Temperli P, Pollo C, Vingerhoets F, Villemure JG. Electrophysiological localization of the subthalamic nucleus in parkinsonian patients. Neurosci Lett 2002;325:144–6. Roby RJ. A simplified circuit for stimulus artifact suppression. Electroencephalogr Clin Neurophysiol 1975;39:85–7. Salisbury JI, Sun Y. Assessment of chaotic parameters in nonstationary electrocardiograms by use of empirical mode decomposition. Ann Biomed Eng 2004;32:1348–54. Ville J. Théorie et application de la notion de signal analytique. Câbles 637 et Transmissions 1948;2A:61–74. Wagenaar D, Potter S. Real-time multi-channel stimulus artifact suppression by local curve fitting. J Neurosci Methods 2002;120:113–20. Welter ML, Houeto JL, Bonnet AM, Bejjani PB, Mesnage V, Dormont D, et al. Effects of high-frequency stimulation on subthalamic neuronal activity in parkinsonian patients. Arch Neurol 2004;61:89–96. Wichmann T. A digital averaging method for removal of stimulus artifacts in neurophysiologic experiments. J Neurosci Methods 2000;98:57–62. Wu Z, Huang NE. Ensemble empirical mode decomposition: a noise-assisted data analysis method. Adv Adapt Data Anal 2009;1:1–41. Zhang F, Miller CA, Robinson BK, Abbas PJ, Hu N. Changes across time in spike rate and spike amplitude of auditory nerve fibers stimulated by electric pulse trains. J Assoc Res Otolaryngol 2007;8:356–72.