Physics Letters A 356 (2006) 237–241 www.elsevier.com/locate/pla

Temporal structure of neuronal population oscillations with empirical model decomposition Xiaoli Li Institute of Electrical Engineering, Yanshan University, Qinhuangdao 066004, China Received 22 December 2005; received in revised form 17 March 2006; accepted 21 March 2006 Available online 31 March 2006 Communicated by C.R. Doering

Abstract Frequency analysis of neuronal oscillation is very important for understanding the neural information processing and mechanism of disorder in the brain. This Letter addresses a new method to analyze the neuronal population oscillations with empirical mode decomposition (EMD). Following EMD of neuronal oscillation, a series of intrinsic mode functions (IMFs) are obtained, then Hilbert transform of IMFs can be used to extract the instantaneous time frequency structure of neuronal oscillation. The method is applied to analyze the neuronal oscillation in the hippocampus of epileptic rats in vivo, the results show the neuronal oscillations have different descriptions during the pre-ictal, seizure onset and ictal periods of the epileptic EEG at the different frequency band. This new method is very helpful to provide a view for the temporal structure of neural oscillation. © 2006 Elsevier B.V. All rights reserved. Keywords: Neuronal oscillation; Empirical mode decomposition; Epileptic seizure; Hippocampus

1. Introduction Analysis of population neuronal oscillations is a fundamental issue in neuroscience [1,2]. This is because the neuronal oscillations are directly related to the dynamic property of interaction between the cellular and synaptic mechanisms [3]. Through the analysis of neuronal oscillation, we can understand the mechanism of information processing or disorder in the neural network, such as the organization and size of the neural network. To give insight on the mechanism of brain information processing, Fourier and wavelet transforms are often used to analyze neural activity (e.g. [4–6]). Fourier transforms can perfectly represent a harmonic oscillation signal, and wavelet transforms can present the local information of a non-stationary signal (e.g. [5,6]). It is known that wavelet transforms based frequency distribution of neural activity depends on the selection

E-mail address: [email protected] (X. Li). 0375-9601/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.physleta.2006.03.045

of basis function. Fourier analysis often need some kind of window to eliminate the end effects, however it is not find a guide to choose the window size during the non-stationary process. EMD in conjunction with Hilbert transform can be used to analyze some complicate signals from biology (e.g. [7–9]). One of the advantages of this new method is that it can break down a complicated signal without a basis function such as sine or wavelet functions. All of the oscillatory modes embedded in the complicated signal can be described in the time frequency domain. The amplitude and time of each oscillatory mode are not distorted due to the basic function, for example, the phase shift. The extracted modes have well behaved Hilbert transforms, thus the instantaneous frequency and amplitude of each mode can be accurately calculated, so a new frequency distribution of EEG data may be obtained. In this Letter, we addresses this new decomposition method to describe the details of neural oscillations in the rat hippocampus during epileptic seizures arising from a focus induced by a single intra-hippocampal injection of a low dose of tetanus toxin [10,11].

238

X. Li / Physics Letters A 356 (2006) 237–241

2. Methods Once an oscillatory function satisfies the following two conditions: the number of extrema and the number of zero either be equal or differ at most by one and the mean value of the envelop of the local maxima and the envelop of the local minima must be zero at any time, it is called as an intrinsic mode function (IMF) [7]. For a given non-stationary signal s(t), the EMD method can decompose the signal as a linear combination of intrinsic mode functions (IMFs), Cn (n = 1, 2, . . . , N ), N is the number of IMFs. The EMD algorithm [7] involves the following steps: Step 1: Set n = 0 and s0 (t) = s(t). Define hk = sn (t), k = 0; Step 2: Identification of all the extrema (maxima and minima) of the series hk (t); Step 3: Generation of the upper and lower envelope, U (t) and L(t), via cubic spline interpolation for all of the maxima and minima, respectively; Step 4: Shift. The mean of U (t) and L(t) is denoted by mk (t) (i.e. mk (t) = (U (t) + L(t))/2), and subtraction of mk (t) to obtain an IMF candidate hk+1 = hk − mk (t). The functions U (t) and L(t) should envelop the data between them, i.e. L(t)  hk (t)  U (t) for all t . If hk+1 is not an intrinsic mode functions (IMF), then increment k, return to Step 2 and take hk+1 (t) in place of

hk (t); else, define the IMF Cn (t) as hk+1 (t) and the residual sn+1 (t) = sn (t) − Cn (t). The IMF should satisfy two conditions: the local maxima of the data are positive, and the local minima are negative and the mean value of the envelopes defined by the local maxima and by the local minima need to be zero for all of the analyzed EEG data. If the convergence criteria are not met, increment n and return to Step 1. Typical convergence criteria are to test whether the residual is either smaller than a predetermined value or a monotonic function. At the end of first shift process, we will gain the highest frequency oscillation and the oscillation is extracted from original oscillation. Next shift process will extract the secondary fast oscillation. By repeating the shift process, it will extract the oscillations embedded with different frequencies and a trend data from the original EEG recording. The instantaneous phase of an oscillatory time series with narrowband can be obtained by means of the analytic signal concept proposed by Gabor in 1946. Given a time series x(t), the complex process of this time series can be written as (t) is x(t) + ixH (t) = A(t) exp(iϕ(t)), where the function  ∞ xH x(t) the Hilbert transforms of x(t), xH (t) = π −1 P.V. −∞ t−τ dτ , where P.V. means that the integral is taken in the sense of the Cauchy principal value. The instantaneous amplitude 1 dϕ (A(t)) and phase ϕ(t) or frequency (f (t) = 2π dt ) of each IMF can be gained by above Hilbert transforms. For all of IMFs Cn (t) (n = 1, 2, . . . , N ) in EEG data, a series of instantaneous amplitudes and frequencies can be obtained, ai (t) and fi (t) (i = 1, . . . , N ). These instantaneous amplitudes and

Fig. 1. (Top trace) the EEG recordings of CA1 in the right hippocampus of epileptic rat from normal to seizure states (duration of this recording is 15 seconds, the sampling frequency is 2.5 kHz). The traces can be decomposed into pre-ictal and ictal. The segments of I and II in the right CA1 are extracted for further analysis. (Right traces) the empirical mode decomposition of an EEG data during the pre-ictal stage. EEG recordings are 2 seconds. The EMD of the EEG data, including 9 components. The components are listed from high to low frequencies. (Left traces) the empirical mode decomposition of an EEG data during the ictal stage. The EEG data is a typical rhythmic oscillation. Fast oscillations are superimposed on slow oscillations. EMD of the EEG data, including 9 components. The amplitude of fast ripples is very small, and they are superimposed to a slow oscillations. The fast oscillations are distributed at the top of a slow oscillation (row 8 around 2 Hz), then a series of fast oscillations (row 4–7) are superimposed on the negative phase of a slow oscillation (row 8).

X. Li / Physics Letters A 356 (2006) 237–241

239

frequencies can form a frequency distribution of the EEG data. 3. Results 3.1. Hippocampal EEG recording In this study, the rat tetanus toxin model of focal epilepsy is applied to study the neuronal oscillations in CA1/CA3 of rat hippocampus. Recordings were made during previous studies of this model [10,11]. Methods were described in the previous reports, but briefly, male Sprague–Dawley rats were anesthetised with halothane, electrodes were placed into CA1 and CA3 of both hippocampi, and 12 mouse LD50 of tetanus toxin (Welcome Foundation Research Laboratories, Beckenham, Kent, UK) was injected in 1 µl saline into the right hippocampus (3.5 mm lateral and 2.7 mm posterior to bregma) over 1 min. Then animals were allowed to recover with free access to food and water. Recording started before the onset of spontaneous epileptiform discharges. The amplified EEG signal was stored on FM tape (Racal V Store, Racal, Southampton, UK, band-pass DC-3.25 kHz). The data were sampled at 2.5 kHz using Spike2 (CED Ltd., UK) running on a personal computer for further analysis. More detail of the EEG recording can be found in [10,11]. In this study, the 9 seizures of 6 rats are investigated, namely, including 36 recordings (left CA1/CA3 and right CA1/CA3 in hippocampal areas, CA1 and CA3 are two main areas of hippocampus, respectively). 3.2. Empirical mode decomposition to EEG recordings The EEG recording in Fig. 1 describes a typical process of pre-ictal towards a seizure. To identity the difference of neuronal oscillations during different stages, we often extract three segments: pre-ictal (I) and ictal (II) from the EEG recordings (see Fig. 1 top trace, R. CA1) for further analysis. Observing Fig. 1, it is found that the neuronal population oscillations are non-stationary, some high frequency signals are associated with a low frequency oscillation. Fig. 1 (right traces) shows the EMD of the pre-ictal activity. The data is decomposed into 9 IMFs (or suboscillations); the frequency range for channels 1–9 is from 200 to 600 Hz (channel 1), 80–200 Hz (channel 2), 30 to 80 Hz (channels 3 and 4), 14– 30 Hz (channel 5), 4–14 Hz (channels 6 and 7), and 0–4 Hz (channel 8 and 9). Fig. 1 (left traces) plots the ictal EEG recording and its EMD. It is clear that the ictal EEG is composed of large regular slow oscillations and regular fast oscillations. The slow oscillations modulate the fast oscillations. The EEG data s(t) can be exactly reconstructed by a linear superposition: s(t) = N i=1 Ci (t) + rN+1 (t), where rN+1 (t) is the residual of the EEG data that presents a trend. EMD method should be a direct extraction of instantaneous amplitude associated with various intrinsic time-frequencies. Following EMD, the instantaneous frequency of IMFs with Hilbert transform (HT) can be calculated; as a result, a time frequency distribution of EEG recording can be obtained. In Fig. 2, the time frequency distribution of two segments in Fig. 1 is plot-

Fig. 2. The spectrum of EEG recording with EMD and Hilbert transform. (A) The time frequency distribution at the pre-ictal stage; (B) The time frequency at the ictal stage; (C) The marginal frequency distribution of two neural activities at the different stages.

ted. In the pre-ictal (seeing Fig. 2A), the frequency distribution concentrate on the low band, the information of high frequency band is very sparse. However, the frequency distribution moves to high frequency band (see Fig. 2B). The marginal spectrums of both EMD-HT are plotted in Fig. 2C, the difference of both is very significant with this new method. To identify the difference of pre-ictal and ictal EEG data with above new method, the instantaneous amplitude and frequency of each IMF of the length of 4 seconds are calculated. Then, the median frequency and its amplitude of each IMF are sorted out. The frequency distribution of an EEG based on EMD can be obtained. Fig. 3 (left) shows the frequency distribution of right CA1/CA3 and left CA1/CA3 based on the 9 seizures from 6 rats. We can make two significant observations. One is that the activity of neurons in CA1 is more significant by comparison with the distribution of frequency at CA1 and at CA3, since the amplitude of neural activity at CA1 is higher and the distribution of frequency is wider. Another is the significant difference at the frequency distribution between the pre-ictal and ictal stages, in particular, at the theta, beta and gamma bands, as can be seen in Fig. 3 (right).

240

X. Li / Physics Letters A 356 (2006) 237–241

Fig. 3. The frequency distribution at the pre-ictal and ictal stage with the empirical mode decomposition and Hilbert transforms. Left figures: x axis is the median frequency of each IMFs, y axis is the amplitude corresponding to the median frequency. Right figures: I: 1–4 Hz, II: 4–12 Hz, III: 12–30 Hz, IV: 30–100 Hz, V: 100–200 Hz, VI: 200–600 Hz; y axis is the mean and standard deviation of amplitudes at each frequency band.

4. Discussions Investigating the long-range interactions involved in seizure generation depends on intact, preferably freely-moving, animal models in vivo. In this study, the tetanus toxin model of epilepsy is applied on the behaving rats. EMD is used to decompose the neuronal oscillations from right and left CA1 and CA3. According to the examination of 9 seizures from 6 rates, it is found that the frequency distribution of neuronal oscillations shows that the difference of dynamic activity between the pre-ictal and ictal stage is remarkable at the beta and gamma waves.

The analysis of oscillations of neural network becomes a very helpful approach to understand fundamental issues in cognitive processing and neural clinics [1,2]. Population oscillations of neural networks are very complicated due to the chemical and electrical interactions (connections) of thousands of neurons in vivo [12]. The analysis of population oscillations is still a challenging topic in neural signal processing. The core in neuronal oscillation analysis is to discover how slow oscillation modulated fast oscillation, which will be useful to indicate the mechanism of neural signal processing. In this Letter, EMD and Hilbert transform is addressed to analyze neuronal population

X. Li / Physics Letters A 356 (2006) 237–241

oscillation. Although the details of neural population oscillation is not very clear until now, the decomposition of neural population oscillation with EMD method seems to be able to indicate a mixture relation of neural oscillations of different frequency. The IMFs with high frequency information correspond to fast neural oscillation; ones with low frequency information correspond to slow neural oscillation. Herein, some interesting oscillations gamma wave (30–100 Hz) and beta waves (12– 30 Hz) are discussed. An attractive idea is that the pyramidal cells and interneurons are phase-locked to generate the gamma wave [13]. As can be seen in Fig. 3, the amplitude gamma wave increase sharply from the pre-ictal to ictal stages. This finding is in agreement with above idea. Some evidences indicate the long rang synchronization of neural networks is modulated by beta oscillation [14]. In Fig. 3, it is found that the beta oscillation change significantly form pre-ictal to ictal stages, the beta wave is derived from the spread of seizure. These findings strongly support the long rang synchronization mechanism to generate beta wave. On the other hand, it is possible that the beta oscillation plays a very important role in the generation procedure of seizure. Furthermore, the significant change of beta wave at CA1 region than at CA3 region indicates the beta generator is related to the connectivity of neural networks, since the connectivity of neural network at CA1 and CA3 is different. Based on the EMD method, the EEG recordings in the hippocampus of the behaving rats are composed of several relaxation oscillations, which supports the idea that the brain dynamics can be modelled as relaxation oscillators, for example, Fitzhugh–Nagumo model (e.g. [15,16]). In the Fitzhugh– Nagumo model, the neurons or neural networks contain their natural frequency and some driving frequency from other neurons or neural networks. When a neuron receives enough electrical input from other neurons it will fire (sending off its own electrical pulse); then the neuron starts to recover by receiving the ions from other neurons). In this study, it is found the behaviour of the relaxation oscillator from EMD decomposition, which is in line with our observations (in Fig. 1). In this Letter, we find that the EMD method can discover the complex relation among several oscillations from the examples of epileptic EEG data. The new method can provide a direct way, without any transform to break down a complicated EEG into several oscillations. From this viewpoint, EMD is superior in analyzing non-linear and non-stationary EEG signals. The frequency distribution of each oscillation based on the EMD is similar to the result with the discrete wavelet transform. However, there is a clear difference between EMD and discrete wavelet transform, EMD is a local and data-driven technique, the basic function are directly derived from EEG data itself [17], unlike discrete wavelet transform it needs to select a filter function. Moreover, it is noteworthy that there also exist some drawbacks on EMD method. The shifting process is based on the time series itself, so the length of time series has an im-

241

portant effect on its decomposition. The aim of this Letter is to indicate the difference of frequency distribution of neuronal population oscillations; the segments analyzed are enough long to avoid the negative effect. However, it is noteworthy that the operation of EMD is similar to one of discrete wavelet transform. IMFs need to satisfy the two conditions, so this method is not very good to be applied to treat with a time series with transient of very short time; IMFs are practically orthogonal within a sufficient diminutive interval at a certain instant in time [18], mathematical orthogonality of the set of IMFs is not guaranteed, since these intrinsic functions are not given in a closed analytical form [19]. In practice EMD just is a pre-process of Hilbert transform, the real physical meaning of time series is derived from the instantaneous frequency based on Hilbert transform. In summary, an EEG data can be well decomposed into a few suboscillations using EMD without time shift and basis function. A frequency distribution of EEG data can be obtained. EMD method can discover the more details of an EEG recording by comparing with some of the traditional methods. Acknowledgements This Letter is supported by National Natural Science Foundation of China (No. 60575012). Author appreciates John Jefferys and John Fox for the EEG data and helpful comments. References [1] G. Buzsaki, A. Draguhn, Science 304 (2004) 1926. [2] A.K. Engel, P. Fries, W. Singer, Nat. Rev. Neurosci. 2 (10) (2001) 704. [3] X.-J. Wang, Neural oscillations, in: Encyclopedia of Cognitive Science, MacMillan Reference, 2003, pp. 272–280. [4] G. Lantz, C.M. Michel, M. Seeck, O. Blanke, T. Landis, I. Rosen, Clin. Neurophysiol. 110 (1) (1999) 176. [5] H. Adeli, Z. Zhou, N. Dadmehr, J. Neurosci. Methods 123 (1) (2003) 69. [6] Y.U. Khan, J. Gotman, Clin. Neurophysiol. 114 (5) (2003) 898. [7] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, Proc. R. Soc. London 454 (1998) 903. [8] R. Balocchi, D. Menicucci, E. Santarcangelo, L. Sebastiani, A. Gemignani, B. Ghelarducci, M. Varanini, Chaos Solitons Fractals 20 (1) (2004) 171. [9] H. Liang, Z. Lin, R.W. McCallum, Med. Biol. Eng. Comput. 38 (2000) 35. [10] G.T. Finnerty, J.G.R. Jefferys, J. Neurophysiol. 83 (4) (2000) 2217. [11] G.T. Finnerty, J.G.R. Jefferys, J. Neurophysiol. 88 (6) (2002) 2919. [12] R. Llinas, Science 242 (4886) (1988) 1654. [13] J. Csicsvari, B. Jamieson, K.D. Wise, G. Buzsáki, Neuron 37 (2) (2003) 311. [14] A. Schnitzler, J. Gross, Neuroscience 6 (2005) 286. [15] R. Fitzhugh, Biophys. J. 1 (1961) 445. [16] J. Nagumo, S. Arimoto, S. Yoshizawa, Proc. IRL 50 (1960) 2061. [17] H. Liang, S.L. Bressler, E.A. Buffalo, R. Desimone, P. Fries, Biol. Cybernet. 92 (6) (2005) 380. [18] W. Huang, Z. Shen, N.E. Huang, Y.C. Fung, Proc. Natl. Acad. Sci. 95 (9) (1998) 4816. [19] M. Dätig, T. Schlurmann, Ocean Engrg. 31 (14–15) (2004) 1783.