Epileptic Seizures Detection Based on Empirical Mode Decomposition of EEG Signals

1. Introduction

Epilepsy is a chronic neurological disorder that affects more than 50 million people world wide, characterized by recurrent seizures (World Health Organization [WHO], 2006). An epileptic seizure is a transient occurrence of signs and/or symptoms due to abnormal excessive or synchronous neuronal activity in the brain (Fisher et al., 2005 & Berg et al., 2010). This electrical hyperactivity can have its source in different parts of the brain and produces physical symptoms such as short periods of inattention and loose of memory, a sensory hallucination, or a whole-body convulsion. The frequency of these events can vary from one in a year to several in a day. The majority of the patients suffer from unpredictable, persistent and frequent seizures which limit the independence of an individual, increase the risk of serious injury and mobility, and result in both social isolation and economic hardship (Friedman & Gilliam, 2010). In addition, the patients with epilepsy have an increased mortality risk of approximately 2 to 3 times that of the general population (Ficker, 2000).

The first line of treatment for epilepsy is with multiple anti-epileptic drugs and it is effective in about 70% of the cases. From the 30% remaining affected individuals only less than 10% could benefit from surgical therapy leaving a 20% of the total of people with epilepsy who will continue suffering sudden, incontrollable seizures and for whom other forms of treatment are being investigated (Theodore & Ficker, 2004; WHO, 2006).

For any of the reasons exposed before the seizure detection is an important component in the diagnosis of epilepsy and for the seizures control. In the clinical practice this detection basically involves visual scanning of Electroencephalogram (EEG) long recordings by the physicians in order to detect and classify the seizure activity present in the EEG signal. Usually these are multichannel records of 24 to 72 hours length which implies a very time consuming task and it is also kwon that the conclusions are very subjective so disagreement between physicians are not rare.

The seek here is to detect automatically in long term EEG records those segments of the signal that present epileptic seizures for the shake of reducing the high amount of information to be analyzed by the neurologists. Thus them could focus their attention in these part of the information so a more precisely and quick diagnosis can be made. Seizure detection is also a useful tool for treatments such us timely drug delivery, electrical stimulation and seizure alert systems.

Automated seizure detection, quantification and recognition have been of interest of the biomedical community researchers since the 1970s. In some initial works a number of parameters of EEG waves such us amplitude, sharpness and duration were measured and evaluated (Sanei & Chambers, 2007). This first approach is sensitive to artefacts so in the following years numerous and diverse techniques have been employed and refined to improved epileptic seizures detection.

Artificial Neural Networks (ANN) have been used both to detect abnormal patterns in the EEG (Schuyler et al., 2007 & Baoet al., 2009) and as seizure parameters classifier (Tzallaset al., 2007). Wavelet Transform is also widely used for epilepsy detection (Adeliet al., 2007). Others studies combine Approximate Entropy and Lempel-Ziv Complexity (Abásolo et al., 2007 & Zandiet al., 2009), and Time Frequency Distributions (Tzallas et al., 2007).

In the studies referenced in the previous paragraphs, it had been proposed different seizure detectors that had been tested in particular EEG databases each. In some cases the epileptic EEG records were of a few seconds long. Other techniques were implemented in rats’ EEGs with induced seizures were also used. The most recently works used long term epileptic EEGs for a small number of patients or grouped by the type of epilepsy they suffer. In this sense due to seizure detection algorithms were not evaluated on the same database to date so no standardization exists about the good performance of an epileptic seizure detector (Varsavsky et al., 2011).

The aim of this chapter is to examine the recent Empirical Mode Decomposition (EMD) technique for the extraction of features of epileptic EEG records to be used in two seizure detectors. The algorithms will be tested in 21 multichannel EEG recordings of patients suffering different focal epilepsies. Along the sections of this chapter it will be described the used EEG records, the EMD algorithm as well as the features extracted to be used in the developed seizures detectors, the obtained results and finally the conclusions and discussion will be exposed.

2. The EEG database

The EEG database contains invasive EEG recordings of 21 patients suffering from medically intractable focal epilepsy. The data were recorded during invasive pre-surgical epilepsy monitoring at the Epilepsy Center of the University Hospital of Freiburg, Germany (Freiburg, 2008). In order to obtain a high signal-to-noise ratio, fewer artifacts, and to record directly from focal areas, intracranial grid-, strip-, and depth-electrodes were used. The EEG data were acquired using a Neurofile NT digital video EEG system with 128 channels, 256 Hz sampling rate, and a 16 bits A/D converter. Notch or band pass filters have not been applied in the acquisition stage.

The available EEG records include only 6 channels (3 focal electrodes and 3 extrafocal electrodes). The records are divided into segments of 1 hour long. In this study, only the 3 focal channels were used. A total of 87 seizures from 21 patients (8M, 13F, age: 29.9 ± 11.9 years) were analyzed. The details of the database are summarized in Table 1.

3.Empirical Mode Decomposition

In the last years, a technique called Empirical Mode Decomposition (EMD) has been proposed for the analysis of non-linear and non-stationary series (Huang et al., 1998). The EMD adaptively decomposes a signal into oscillating components or Intrinsic Mode Functions (IMFs). The EMD is in fact a type of filter bank decomposition method whose sub bands are built as needed to separate the different natural components of the signal. In the field of biomedical signal processing EMD has been used for the analysis of respiratory mechanomyographic signals (Torres et al., 2007), for denoising in ECG records (Benget al., 2006). Particularly, this technique was implemented to extract features from EEG signals for mental task classification (Diezet al., 2009), it was used to obtain adaptive bands on EEG signals (Diezet al., 2011) and also for epileptic seizure detection in EEG signals in 5 patients with temporal lobe focal epilepsy (Tafreshiet al., 2008). In this sense the authors of this chapter have previously developed algorithms based on EMD for seizure detection and they have been tested in 9 long EEG records of patients with temporal focal epilepsy (Oroscoet al., 2009) and in 21 patients with different epilepsies (Oroscoet al., 2010).

Table 1.

Freiburg EEG Database.

3.1. The EMD algorithm

The EMD is a general nonlinear non-stationary signal decomposition method. The aim of the EMD is to decompose the signal into a sum of Intrinsic Mode Functions (IMFs). An IMF is defined as a function that satisfies two conditions (Huang et al., 1998):

1. In the entire signal, the number of extrema and the number of zero crossings must be equal or differ at most by one.

2. At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima must be zero (or close to zero).

The major advantage of the EMD is that the IMFs are derived directly from the signal itself and does not require any a priori known basis. Hence the analysis is adaptive, in contrast to Fourier or Wavelet Transform, where the signal is decomposed in a linear combination of predefined basis functions.

Given a signal x(t) such us it is showed in figure 1, the algorithm of the EMD can be summarized in the following 6 steps (Huang et al., 1998):

1. Find local maxima and minima of d₀(t)=x(t).

2. Interpolate between the maxima and minima in order to obtain the upper and lower envelopes e_u(t) and e_l(t), respectively.

3. Compute the mean of the envelopes m(t)=(e_u(t)+ e_l(t))/2.

4. Extract the detail d₁(t)= d₀(t)-m(t)

5. Iterate steps 1-4 on the residual until the detail signal d_k(t) can be considered an IMF (accomplish the two conditions): c₁(t)= d_k(t)

6. Iterate steps 1-5 on the residual r_n(t)=x(t)- c_n(t) in order to obtain all the IMFs c₁(t),.., c_N(t) of the signal.

The procedure terminates when the residual c_N(t) is either a constant, a monotonic slope, or a function with only one extrema.

The result of the EMD process produces N IMFs (c₁(t), …, c_N(t)) and a residue signal (r_N(t)):

$$x(t)={\displaystyle \sum _{n=1}^N{c}_n(t)+r_N(t)}$$E1

Figure 2 shows the complete process of EMD for the example signal x(t). It can be observed that the lower order IMFs capture fast oscillation modes of the signal, while the higher order IMFs capture the slow oscillation modes.

The EMD is a technique essentially defined by an algorithm and there is not an analytical formulation to obtain the IMFs. Furthermore, several algorithmic variations have been proposed in order to obtain the IMFs decomposition. In this work it had been used the algorithm proposed by Flandrin (2007) & Rillinget al. (2009), in which, in order to accomplish the second IMF condition, it is utilized a criterion that compares the amplitude of the mean of the upper and lower envelopes with the amplitude of the corresponding IMF. This criterion is based on two thresholds (θ₁ and θ₂) and a tolerance parameter (α). It were also used the default values proposed by Rillinget al. (2009): α=0.05, θ₁=0.05 and θ₂=0.5.

3.2. EMD applied to EEG analysis

For the purposes of this work the EMD of the EEG signals was achieved computing IMF1 to IMF5 for every segments of each channel. After several initial tests it was concluded that IMF4 and IMF5 do not contributed to seizure detection, so they were discarded. Thus IMF1, IMF2 and IMF3 of each segment of EEG signals were used in further analysis.

Figure 3 shows an example of a 300 s EEG segment without seizure for one channel and their first 3 IMFs obtained with the described EMD method. Figure 4 illustrates a 300 s EEG segment with an epileptic seizure of the same patient and their corresponding first 3 IMFs. In figure 3 it can be observed how the energy of the IMF remains approximately between the same levels along the showed time period while for the EEG segment of figure 4 the mode functions highlight the increased energy during the seizure.

Figure 1.

Step 1, 2, 3 and 4 of the EMD algorithm. In the top pannel the original signal, in the middle pannel the upper (blue) and the lower (red) envelopes are showed as well as the mean of them (magenta). In the bottom pannel the obtained residue. The figure is a modified reproduction of figures available in http://perso.ens-lyon.fr/patrick.flandrin/emd.html

Figure 2.

Original signal x(t) and the result of its EMD computation. The figure is a modified reproduction of figures available in http://perso.ens-lyon.fr/patrick.flandrin/emd.html

Figure 3.

EEG segment without seizure for one channel and IMF1 to IMF3 of the signal.

Figure 4.

EEG segment with a seizure for one channel and IMF1 to IMF3 of the signal. Red lines indiate the seizure time endpoints established by the neurologists.

4. Features and detectors

In this chapter two different epileptic seizure detectors based on the EMD of EEG signals will be described. In the first detector, the algorithm computes the energy of each IMF and performs the detection based on an energy threshold and a minimum seizure duration decision. The second detector consists on the extraction of several time and frequency features of IMFs, subsequently a feature selection based on a Mann-Whitney test and Lambda of Wilks criterion is performed and in a last stage linear discriminant analysis (LDA) of the selected parameters is used to classify epileptic seizure and normal EEG segments. In figure 5 the block diagrams of both detectors are showed.

Figure 5.

Block diagrams of two epileptic seizure detectors.

4.1. Preprocessing and EMD

All EEG records were initially filtered with a second order, bidirectional, Butterworth, 50 Hz notch filter in order to remove the power line interference. Then, the EEG signals were band-pass filtered with a second order, bidirectional, Butterworth filter with a bandwidth of 0.5 - 60 Hz.

Next, all EEG records were resampled to 128 Hz in order to reduce computation time of EMD decomposition. This operation does not have any influence on the results since the bandwidth of the signal of interest does not exceed the 60 Hz.

Finally, the EMD of EEG signals was achieved as described in section 3.2.

4.2. First detector

The detector presented in this section and schematized in the left side of figure 5 can be separate in 4 main blocks. The first and second stages consist on the preprocessing of EEG signal and the EMD computation as described in 4.1. The third stage implies the energy computation and the last one, and the most complex, is the seizure detection strategy itself.

4.2.1. Energy computation

The first proposed algorithm takes the IMF1, IMF2 and IMF3 of the EEG signals of each channel and computes the energy serie (ENi) of each IMFi as shown in (2).

$$ENi(n)=\frac{1}{L}\sqrt{{\displaystyle \sum _{m=n-L/2}^{n+L/2-1}{(IMFi(m))}^2}}i=\mathrm{1,}\mathrm{2,}3$$E2

In equation (2) idenotes the i-th IMF, n is its sample number and L is the length in samples for the energy computing window. In this work a 15 s moving, overlapped window (L=1920 samples) is used. Thus, once this computation ends three energy series (EN1, EN2 and EN3) for all EEG segments of each channel are obtained (see Figures 6 and 7).

4.2.2. Seizure detection method

In first place it will be describe what is called as an event. An event is define here like the energy series portions that overcomes a certain threshold for more than 30 s. The threshold is computed as (3)

Thr_ENi = mean (ENi) + 1.5*std (ENi)

where mean(ENi) and std(ENi) are the mean and the standard deviation values of the i-th energy serie considering the whole EEG channel.

Thus the first stage in this seizure detector is determined all the events present in each energy series of each channel.

The second decision step is identifying those events present in at least two of the three ENs of each channel. This criterion is used in order to discard possible artifacts that could be present in only one ENi.

Finally, in a third stage an interchannel decision is done by choosing the events (selected in the previous stage for each channel) that are present in at least two of the three studied channels.

Hence all events that satisfy the three decision stages are detected as epileptic seizures.

Figure 6.

A no seizure EEG segment of one channel and EN1 to EN3 series of the signal EMD

Figure 7.

A seizure EEG segment of one channel and EN1 to EN3 series of the signal EMD

Figure 6 illustrates the energy series (ENi) of the IMFs showed in figure 3. It is observed that the energy for each of the three IMFs does not overcome its corresponding threshold (computed by equation (3)), which is indicated with the red dashed line. So no event is detected in this EEG channel as well it is not detected in neither of the other channels (that are not showed here) so no seizure is present for this segment, corresponding with the database information.

In figure 7 the energy series (ENi) of the IMFs showed in figure 4 are illustrated. In this case the energy rise above the threshold in the 3 IMFs and lasts more than 30 s satisfying thus the event condition for each case so for this channel the second decision step is also accomplished. If this occurs for at least one of the remaining channels then a seizure is detected. In this example the events are detected in the three channels and also match with the seizure time endpoints established by the neurologists.

4.3. Second detector

In the right side of figure 5 a block diagram of the second detector is illustrated. In this case the preprocessing stage and the EMD computation (describe in Section 4.1) are the same as the first detector.

Next several time and frequency features of the IMFs are computed and then selected using a Mann-Whitney test and Lambda of Wilks criterion. Finally, a linear discriminant analysis (LDA) is performed to discriminate epileptic seizures and normal EEG segments.

4.3.1. Feature extraction

In order to characterize the EEG signals several features were computed upon these 3 IMFs series (IMF1 to 3) calculated for each channel. For each IMF, a set of parameters in time and frequency domains were computed.

In this stage in order to improve the statistical stationary of EEG records each IMF was divided in segments of 15 s. Hence the whole IMFs selected of the all EEG records analyzed computes a total of 45517 segments, 44828 of them without epileptic seizures and 689 segments denoted as having only one epileptic seizure each.

In time domain, the following parameters were calculated on each IMF: coefficient of variation (VC), Median Absolute Deviation (MAD), Standard Deviation (STD), Mean Value (MV), Variance (VAR) and Root Mean Square Value (RMS). They are summarized in table 2.

For frequency domain, the power spectral density (PSD) of IMF1, IMF2 and IMF3 was estimated by the periodogram method with a Hanning window.

Then, classical parameters of descriptive statistics were computed on the PSD. Therefore, the following frequency features were obtained on the spectrum of each IMF: Central, Mean and Peak Frequencies (CF, MF and PF), Standard Deviation Frequency (STDF), First and Third Quartile Frequencies (Q₁F, Q₃F), Interquartile Range (IR), 95% cumulated energy Frequency (MAXF), Asymmetry Coefficient (AC) and Kurtosis Coefficient (KC) (Marple, 1987). These frequency parameters are listed in table 3.

Table 2.

Time Domain Features

Table 3.

Frequency Domain Features

Resuming, 10 frequency domain parameters and 6 time domain features were computed. Thus, for IMF1 we have 16 parameters for each 15 second segment obtaining in this way a series with the time evolution of each feature. The same procedure is repeated for IMF2 and IMF3. Hence this implies a computation of 48 features series for each EEG channel and a total of 144 series considering the three EEG channels.

4.3.2. Feature selection

In order to reduce the dimensionality problem, the median of the individual values of each features series for the three channels were initially computed. For example, we take CF of IMF1 of channel 1, CF of IMF1 of channel 2 and CF of IMF1 of channel 3 and calculate the median of this parameter resulting in one series for this feature in IMF1. The procedure is repeated for all the parameters and IMFs. Thus, the number of the total features series is reduced to 48.

Even though the vector of features was reduced, its dimension is still too large. As a second approach, a stepwise method based on the statistical parameter Lambda of Wilks (WL) is performed. In an n-dimensional space constructed with n variables and with the matrixes B_nxn and W_nxn representing the square sum and cross products between groups and within-groups, respectively; the WL can be defined as the ratio between their determinants (Tinsley & Brown, 2000) as it can be see in :

$$WL=\frac{\left|W\right|}{\left|W+B\right|}$$E3

In other words, the WL measures the ratio between within-group variability and total variability, and it is a direct measure of the importance of the variables. Therefore, the most important features for the analysis should be selected, i.e. the variables (features) that contribute with more information. Besides, the correlated variables are discarded in this process (Tinsley & Brown, 2000).

With the aim of contrasting significant differences between groups, the value of WL is transformed into the general multivariate statistical F. If F value for a variable is higher than 3.84 (F to get in) this is included in the analysis and once accepted the variable is rejected if its F value is smaller than 2.71 (F to get out).

Once the WL criterion was applied the features selected were 11, their mean and standard deviation values are summarized in Table 4.

4.3.4. Classification

To detect the EEG segments with epileptic seizure a linear discriminant analysis (LDA) was implemented using the classification functions h. These functions are a linear combination of the discriminant variables (X_m) which allows maximize the differences between groups and minimize the differences within-group and are calculated as (Gil Flores et al., 2001):

$$h_k(q)=b_{k0}+b_{k1}{X}_1(q)+\dots +b_{km}{X}_m(q)$$E4

wherek represents the classification groups, i.e., for seizure and no seizure classes (k =2), m is the quantity of features (in this work, m =11) and q is the case to classify. The computing of b coefficients is showed in equations (6) and (7) (Gil Flores et al., 2001).

$$b_{ki}=(n-g){\displaystyle \sum _{j=1}^q{a}_{ij}}{\overline{X}}_j{}_{n=samplesize}^{g=quantityofgroups}$$E5

$$b_{k0}=-0.5{\displaystyle \sum _{j=1}^q{b}_{kj}}{\overline{X}}_j$$E6

Table 4.

Selected features

For LDA the 50% of data was used as training group and the rest as validation group. Then, a second test was done inverting the training and validation groups. The results are exposed in Section 6.2 using the mean value of SEN and SPE obtained in the validation phase for the two classification tests.

Let g₁ be the seizure group and g₂ the no seizure group, once the classification functions were computing for each group the classification is done satisfying the following criteria:

If h₂(q) > h₁(q) then case q belongs to g₂ otherwise if h₂(q) < h₁(q) case q belongs to g₁.

5. Results

In this section it will be expose the performance of both proposed seizures detectors. In order to evaluate the achievement of the algorithms the following diagnostic categories were considered on the detection stage: true negative (TN), false positive (FP), true positive (TP), false negative (FN). The obtained values for these indexes are contrasted with the segments indicated in the database as having seizure or no seizure by the neurologists. Then the statistical diagnostic indexes of sensitivity (SEN) and specificity (SPE) were also computing (Altman, 1993). These indexes are defined as follows and stated in equations (8) and (9).

Sensitivity (SEN): Is the proportion of epileptic seizures segments correctly detected by the algorithm.

$$SEN(\%)=\frac{TP}{TP+FN}\times 100$$E7

Specificity (SPE): Is the proportion of segments without seizures correctly identified by the algorithm.

$$SPE(\%)=\frac{TN}{TN+FP}\times 100$$E8

5.1. Results of first detector

As a first approach for this algorithm its performance was evaluated in two ways. In first place the detector was tested on the data sorted by epilepsy types and then the EEG signals were evaluated all without a specific arrange.

Then in table 5 are resumed the statistical diagnostic indexes of SEN and SPE computed for the different types of epilepsies individually and for the epilepsies all together.

Table 5.

Statistical diagnostic indexes of SEN and SPE for first detector.

5.2. Results of second detector

In order to improve the results obtained with first detector, the second detection scheme detailed in section 4.3 were tested in the same EEG records. Table 6 shows the mean value of SEN and SPE obtained in the validation phase for the two classification tests described in section 4.3.4.

Table 6.

Statistical diagnostic indexes of SEN and SPE for second detector.

6. Discussion and conclusions

Epileptic seizure detection in EEG records is a useful and important tool due to their various applications such us epilepsy research treatments like timely drug delivery, electrical stimulation and seizure alert systems besides diagnostic applications. In this sense it is a real need the development of automatic algorithms that could be able to detect seizures independently of its brain source. It is also important to establish some kind of standardization of the detectors using to test them the same database so a robust comparison of their performance could be carried out.

In this chapter two epileptic seizure detection methods based on the Empirical Mode Decomposition (EMD) of EEG signals has been proposed. On one hand, the use of EMD for seizures detection it is a recent approach. In addition, as a contribution to the setted out problem, long term epileptic EEG intracranial records with different focal epilepsies are used to evaluate the performance of both seizures detectors.

The used EMD algorithm in this work is the one proposed by Flandrin (2007) & Rillinget al. (2009). This technique seems to be more suitable for epileptic EEG records than others of the signal processing area due to the EEG signal presents nonlinear and non-stationary properties during a seizure. Nevertheless, it was recently reported for this version of the algorithm the problem of what is called mode mixing so to solve this a new approach known as Ensemble EMD (EEMD) has been proposed (Wu & Huang, 2009). There are also some extensions of standard EMD to multivariate signals defined by Rehman & Mandic (2010) as Multivariate EMD. Even though the EMD showed a relatively good performance in seizure detection it was observed that the computation time of EMD for each segment is quite time-computing extensive which could represent a disadvantage for analyzing long EEG records. It can be noted that the proposed EMD technique has still much aspects to explore and innovate so its performance could be further improve.

In order to have a complete evaluation of the detectors’ performance they both were first tested making a discrimination of the EEG signals by epilepsy type and then the data were used all without a specific arrange. For the first detector the values of SPE obtained were high, arising up to 90% for the epilepsies grouped like “Others” while the SEN results were non-satisfactory, been 56.8% the highest value for temporal lobe epilepsy records. Whereas the performance of this detector for the complete set of data showed a global SEN and SPE values of 41.4% and 79.3%, respectively.

The results shown in Table 6 indicate that the second detector have remarkably improved the SEN values compared with those obtained for the first detector for all classes of epilepsies. With respect to the ESP values, the results of the second detector were better for temporal lobe epilepsy signals and decrease slightly for the remaining classes of epilepsies. So the global performance of the second detector (SEN = 69.4% and SPE = 69.2%) can be considered satisfactory better than the first one because both values are in the same order.

Other authors had also recently used the Freiburg´s database for seizure detection so a comparison of their works with ours could be made. Henriksenet al. (2010) in their research uses features of Wavelet Transform (WT) of 16 patients (instead of the 21) of the database and classified them by a support vector machine in order to implement an automatic seizure detection algorithm. They obtained a global SEN of 86% and a false detection rate of 0.39/h, but the SPE value is not reported. In a recent work Vardhan & Majumdar (2011) introduce a differential operator to accentuate the seizure part of depth electrode recordings (ECoG) relative to the non-seizure one. The technique was only applied to 5 patients of Freiburg´s database. For 4 patients, they reported 18 of 20 true detections and 2 false detections. The results for the remaining patient were not reported. Finally, Chua et al. (2011) applied the Gotman algorithm (Gotman, 1999) to 15 patients of the database and it was adjusted for epilepsy type with the aim of improve the off-line automated seizure detection methods that will decrease the workload of EEG monitoring units. The obtained values were 78% of SEN and a true positive rate of 51%.

Summarizing, even though some of the detectors described in the previous paragraph obtained higher values of SEN than the ones developed in this chapter it has to be said that all the referenced cases use selected records of the database while the authors of this chapter had tested their algorithms using all 21 EEG recordings available in Freiburg database.

It may also be highlighted that the values of SEN and SPE of first and second detectors could be improved in order to obtain a more reliable application. In this sense, more tests and some adjustments on the algorithms must be made done to be suitable for medical diagnosis. It could be concluded that the developed methods based on EMD are promissory tools for epileptic seizure detection in EEG records.

7. Future works

As future extension of this research in first place the EMD computation time must be reduced may be taking time windows of few seconds to calculate it instead of 1 h EEG segments. It is also needed to improve the values of SEN and SPE so more effort on the features and classifiers must be done.

Acknowledgments

This work has been supported by grants from AgenciaNacional de PromociónCientífica y Tecnológica (ANPCYT - PICT 2006-01689) and Universidad Nacional de San Juan, both institutions from Argentina. The first author is supported by ANPCYT, whereas the second and third authors are supported by CONICET of Argentina.