Fronto-Temporal Analysis of EEG Signals of Patients with Depression: Characterisation, Nonlinear Dynamics and Surrogate Analysis

2.1 Measurement protocol

The real time data was recorded from 30 medication free outpatients under the age group of 20–50 years comprising of 16 female and 14 male patients from the Psychiatry department of Medical College, Calicut, Kerala, India (female mean age: 33; male mean age: 35). The measurement was done on unipolar depression patients who did not have any history of substance abuse and no significant medical illness. Similarly 30 age and sex matched healthy controls also participated in the study who were free of medical illness. None of them reported of a history of any central nervous system disorder. Bipolar EEG recordings using a 24-channel EEG measuring instrument was carried out at locations FP1-T3 (left half) and FP2-T4 (right half) of the brain. The electrodes were placed, based on the International 10–20 electrode placement system. An ear clip electrode attached to the right earlobe served as an isoground connection. The EEG recordings were done by placing the electrodes on the frontopolar-temporal regions both on the left and right half of the brain for duration of 5 minutes each, under eyes closed and eyes open condition in a resting state. The sampling frequency of the signal is 256 Hz and is notch filtered at 50 Hz to remove the power line interference. Statistical analysis was performed by One-way ANOVA to test for differences among the two classes of EEG signals recorded. Informed written consent was obtained from all the subjects who participated in the study and medical ethical committee approval was taken prior to the study. Figures 1 and 2 show a typical EEG signal of normal and depression patient respectively.

2.2 Preprocessing

Artefacts such as eye movements, eye blinks, head movements, cardiac and muscle activation artefacts, tongue movements and power line noise pose a problem for the proper EEG interpretation and analysis. Other artefacts that disrupt the EEG signal include instrument artefacts (faulty electrodes), sweat artefacts, impedance fluctuations, cable movements, pulse artefacts etc. Power line interferences are removed from the EEG signal by using a 50 Hz notch filter. Eye movement and muscle movement artefacts are manually removed from the signal with the help of an expert by visual inspection. In this work, the high frequency components present in the acquired EEG signals are denoised using Total Variation Filtering (TVF) [2].

The TVF employed in this work is based on the algorithm developed by Chambolle [3]. A dual formulation approach is used to minimise the objective function of the Total Variation (TV) denoising problem. So the TV denoising problem amounts to minimising the following discrete function [4]:

where Ais a matrix of size M × N. Smoothing of the signal is controlled by λ, which is known as the regularisation parameter. Since the amount of high frequency noise present in the EEG signal recorded from depression patients using the 24-channel equipment was already low, the optimal value of λfor denoising was found to be 0.9.

2.3 Selection of wavelet

From an array of discrete orthogonal wavelets, Daubechies-1(D1) to D10, Coiflet-1 to Coiflet- 5 and Symlet 1 to Symlet 8, the task is to identify the wavelet which suits well with the individual EEG signal recorded from depression patients. This is necessitated as it is found that there is extreme patient variability and also variability of the signals with respect to the location on the scalp from person to person. All the above 23 wavelets were tested on all the 30 patient records under four categories namely eyes open and eyes closed conditions recorded from the left and right half of the brain.

The best wavelet is chosen based on the highest value of correlation coefficient which indicates a better match of the characteristics of the EEG signal of depression patient with the wavelet selected. Of the 30 cases considered in this experimental study described here, for 85% of the cases, Coiflet 5 emerged as the best suited wavelet. Hence Coiflet 5 is used for analysis.

2.4 Wavelet Transforms

Wavelet Transforms are efficiently used in many of the signal processing applications as it gives more accurate time and frequency representation of the signal. Wavelet Transforms help in the extraction of wavelet coefficients of discrete time signals. An important feature of Wavelet Transform is that it gives good frequency resolution over a large window while good time resolution at high frequencies. This feature has been of great interest to biomedical applications, as most of the biosignals contain high frequency components in a short span and low frequency components over large span. Wavelet Transforms are thus highly useful for the analysis of nonlinear, nonstationary signals as it gives an excellent time–frequency resolution.

In this method, the signal is decomposed into a set of basis functions called wavelets by the Wavelet Transform. These basic functions are obtained by dilations, contractions and shifts of a unique function called wavelet prototype. Continuous wavelets are functions generated from one single function by dilations and translations of a unique mother wavelet ψ (t):

where ais the scale parameter, b,the shifting parameter and t,the time. The function set Ψa,btis called the wavelet family. The Wavelet Transform usually used in engineering application is the Discrete Wavelet Transform (DWT). It uses the discrete values of the scaling and the translational parameters given by,

a=a0jand b=kb0a0jwhere jand kare integers. Then we get:

where jindicates frequency localization and kindicates time localization. Dyadic scheme implementation is the basis for Multiresolution Analysis (MRA) in Discrete Wavelet Transforms. Any time series can be decomposed in terms of coarse approximations provided by scaling functions and the detail information by the wavelet functions [5]. The scaling function is associated with low-pass filters (LPF) and the wavelet function is associated with the high pass filters (HPF). The decomposition of the signal into the different frequency bands is obtained by successive convolution with high-pass and low pass filtering of the time domain signal.

The approximations are the low frequency components and the details are the high frequency components of the time series. The detail coefficients and approximation coefficients at level 1 (CD1 and CA1) are obtained by decimating the outputs from both the filters by 2. The procedure is then repeated by sending the approximation coefficients to the second stage. This is continued till the signal is decomposed at the expected level. In this work, an eight level decomposition is carried out. The EEG signal acquired from the depression patients are sampled at a frequency of 256 Hz. The multiresolution decomposition offers a time-frequency decomposition of the signal involving not only its energy but also the morphological aspects that are relevant for signal recognition and understanding [6]. Each of the wavelet scales corresponds to a specific frequency band given by

where, fis higher frequency limit of the frequency band represented by decomposition level m, f_sis sampling frequency and 2ⁿis the number of data points in the signal. CD1 (64-128 Hz) and CD2 (32-64 Hz) correspond to gamma band, CD3 (16-32 Hz) corresponds to beta band, CD4 (8-16 Hz) corresponds to alpha band, CD5 (4-8 Hz) corresponds to theta band while CD6 (2-4 Hz), CD7 (1-2 Hz), CD8 (0.5-1 Hz) and CA8 (0–0.5 Hz) correspond to delta band. The DC level of the signal corresponding to eighth level approximation coefficient (CA8) is also considered in this work. The signals are reconstructed from the wavelet coefficients for each scale by applying inverse transform. The reconstructed signal coefficients obtained from eight levels of details and from the eighth level of approximation are further used for energy calculations.

2.5 Calculation of energy

After multilevel decomposition, Parseval’s theorem is employed to calculate the energy of the reconstructed signal coefficients at the detail and approximation levels. It helps in identifying the energy distribution in different frequency bands of gamma, beta, alpha, theta and delta [7]. Parseval’s theorem is mathematically expressed as [8]:

where x(n)is the time domain discrete signal, Nis the total number of samples in the signal, ∑n=1Nxn2is the total wave energy of the signal x(n),∑n=1Najn2is the total energy concentrated in the level ‘j’ of the approximated version of the signal. ∑j=1m∑n=1Ndjn2is the total energy concentrated in the detail version of the signal, from level 1 to mand mis the maximum level of wavelet decomposition.

Energy distribution calculations are done on all recordings of eyes open and eyes closed conditions, acquired from both the left (FP1-T3) and right (FP2-T4) frontopolar-temporal regions of the brain for both the normal controls and depression patients. Energy distribution associated with the various bands of frequency, for different levels of detail are plotted in Figures 3–7, for a single case of measurement recorded (from the left half of the brain under eyes closed condition). Similar variations are also observed for measurements with the other recording protocols.

It is observed that there exists a clear difference in the energy levels of EEG signals of both normal and depression patients in the Gamma band (D1) (Figure 3). Values of energy obtained for normal controls are always higher than that of depression patients, covering all age groups. These differences in energy levels tend to narrow down as we move from Gamma band D1 to alpha band D4 through D2 and D3. It is interesting to note that normal subjects register higher energy distribution levels for cases considered up to theta band D5. This trend appears to get reversed beyond theta band D5 (see Figure 7 for Delta band D8). The energy distribution of depression patients in theta band (D5-Figure 6) is almost at par with that of normal controls. But the energy distribution of depression patients has crossed the threshold of normal subjects and is higher than normal in all cases of delta band from D6-D8 and A8. It gives a clear indication that the brain activity of depression patients in gamma, beta and alpha band is lower when compared to healthy controls. Similar is the trend for all other measurements taken from the right side of the brain under eyes closed and open condition and also for measurements from the left half of the brain for eyes open condition.

2.6 Relative Wavelet Energy (RWE)

Relative Wavelet Energy (RWE) gives information about the relative energy associated with the frequency bands and it can detect the degree of similarity between segments of a signal [9, 10]. It is also known from previous studies that RWE is a good tool for detecting and characterising specific phenomenon in time and frequency planes [11].

The energy wavelet coefficients d_j,krepresent the detail signal energy at each resolution level, jgiven by,

The energy scaling coefficients C_kis defined as the energy at decomposition level N + 1given by

Thus the total energy of the signal for all levels is given by

and hence the Relative Wavelet Energy (RWE) is defined as

Clearly ∑jρj=1and the distribution {ρj} can be considered as a time scale density. This provides information to characterise signal energy distribution at different frequency bands.

RWE calculations are carried out on the different frequency bands to understand the variations in healthy subjects and depression patients. Wavelet energy in the gamma band, particularly D1 is negligible, while D2 is negligibly small for depression patients in comparison to normal controls. The RWE levels associated with beta (D3–18.5%) and alpha (D4–22%) bands for normal controls are approximately 25% higher than the corresponding values for depression patients with D3 being 1.9% and D4–4.3% (Figures 8 and 9). RWE of theta band also shows similar trend with normal controls having values slightly lower than twice the values associated with depression patients. RWE associated with delta band (D6–19.5%, D7–23.5%, D8–17.1% and A8–23.7%) show appreciably higher values for depression patients when compared to normal controls. The feature worth noting is that RWE of approximation level has the highest percentage in depression patients to that of normal controls. From the results of RWE calculations, it may be observed that the RWE is more prominent in depression patients in the frequency bands, 1-2 Hz and 0–0.5 Hz. Hence depression may be classified as a very low frequency phenomenon.

RWE values plotted in Figures 8 and 9 represent the calculations carried out on the frequency bands of normal and patient EEG signals recorded from the left half of the brain under eyes open condition. RWE calculations are also carried out on the frequency bands of the EEG signal acquired from the left half of the brain under eyes closed condition and on the EEG signals from the right half of the brain both under eyes open and eyes closed conditions. The observation from all the protocols reveals a high value of RWE in alpha band (D4) of normal controls indicating high activity in the thought process of healthy subjects. Also a high value of RWE in 8th level approximation followed by detail level D7 is observed in depression patients. Hence from all the four cases, it may be concluded that depression phenomenon is confined to the lower frequency bands especially in delta band of 0–4 Hz. In order to analyse the statistical nature of the measurement among the two broad classes of EEG signals recorded, One-way ANOVA is carried out which gave a statistical significant difference (p < 0.005).

2.7 Wavelet Entropy (WE)

The degree of order/disorder associated with a multifrequency signal response is characterised by Wavelet Entropy (WE). The time evolution of WE was also calculated which gave information about the dynamics in the EEG records. It was observed that in contrast to spectral entropy, WE is capable of detecting changes in a nonstationary signal due to the localisation characteristics of the Wavelet Transform. The computational time of WE was significantly shorter since the algorithm involved the use of fast wavelet transform in a multiresolution framework. The results demonstrated that WE could differentiate between specific physiological brain states under spontaneous or stimulus-related conditions. The time evolution of the WE is calculated to give information about the dynamics in the EEG records. A signal generated by a totally random process can be taken as representative of a very disordered behaviour.

The Shannon Wavelet Entropy (WE) as a function of time is calculated [10] as:

where mis the wavelet decomposition level from level 1 to level m.

WE may be considered as a meaningful indicator since it is able to differentiate physiological brain states under normal and depression conditions. Figure 10 represents the Wavelet Entropy calculated for the EEG signals recorded from the left half of the brain under eyes open condition from normal and depression patients. Significant decrease in the WE is observed in the EEG signals of depression patients under all the four recording protocols, indicating a more rhythmic and ordered behaviour of the EEG signal [12, 13]. Being independent of the amplitude or the energy of the signal, the WE gives additional information about EEG signals in comparison to those obtained by using frequency analysis or other standard methods. The use of such quantifiers based on time-frequency methods can contribute to the analysis of brain responses and may also lead to a better understanding of their dynamics.

2.8 Nonlinear measures

Reduction in complexity in patients with disease is the main hypothesis that is checked in most of the research work. Here we analyse the EEG signal complexity and irregularity of the frontopolar-temporal regions of the brain of controls and patients with unipolar depression under resting states of eyes open and closed conditions. The analysis is carried out using Approximate Entropy, Fractal Dimension and Largest Lyapunov exponent.

2.8.1 Approximate Entropy (ApEn)

Approximate Entropy (ApEn) is a statistic quantifying regularity and complexity which has potential application to a wide variety of physiological and clinical time series data. ApEn is a statistical parameter that measures the predictability of the current amplitude values of a physiological signal based on its previous amplitude values. Approximate Entropy is the is the probability difference of the pattern similarities of the connected straight lines of the madjacent points and m + 1 adjacent points of time sequence data [14]. The more complex the sequence data, higher is the probability that new pattern appears and larger the corresponding ApEn. ApEn measures the (logarithmic) likelihood that runs of patterns that are close remain close on next incremental comparisons.

For a time series of Ndata points, {u(i):1 ≤ i ≤ N}, form vector sequences x(1) through x(N-m + 1), defined by x(i) = [u(i),..., u(i + m-1)]. These vectors represent mconsecutive uvalues, commencing with the i^th point. mis the length of compared runs. The distance d[x(i), x(j)] between vectors x(i) and x(j) is defined as the maximum difference in their respective scalar components. Let B_ibe the number of vectors x(j) within rof x(i) for a window length mand let A_ibe the number of vectors x(j) within rof x(i) for a window length m + 1, where ris the tolerance for accepting matches. The function Cimris defined as:

Cimr=BiN−m+1E11

In calculating Cimr, the vector x(i) is called the template and the instance where a vector x(j) is within rof it is called a template match. Cimris the probability that any vector x(j) is within rof x(i). It measures within a tolerance r,(r = k*standard deviation)the regularity, or frequency, of patterns similar to a given pattern of window length m. The function Φmris defined as:

Φmr=1N−m+1∑i=1N−m+1lnCimrE12

where, Φmris the average of the natural logarithms of the functions Cimr.

For finite data sets,

ApEnmrN=Φmr−Φm+1rE13

The parametres N, mand rmust be fixed for each calculation and reffectively works as a filter. The values of mand kadopted in this work are 1 and 0.2 respectively.

The column plot of Figure 11 represents the ApEn values calculated for the EEG signals of healthy controls and patients, acquired from the left half of the brain under eyes open condition. The results of ApEn calculated for from the left part of the brain under eyes closed conditions and from the right part of the brain both under eyes open and closed conditions indicate that normal controls have a higher value of ApEn than depression patients. A low value of ApEn indicates predictability and regularity in a time series, whereas a high value of ApEn indicates unpredictable and random variation. The results of ApEn indicate that the complexity of the brain is high in normal controls while the signals from depression patients are less complex.

2.8.2 Fractal Dimension (FD)

FD analysis is frequently used in biomedical signal processing, including EEG analysis. Fractal Dimension can be used to quantify the complexity and the self-similarity of an object. The EEG FD can expected to be always between 1 and 2 since the dimension of a plane is equal to 2 and the dimension of a line is equal to 1. FD can be calculated using Higuchi’s method, Katz algorithm, box counting approach and so on. FD analysis provides a fast computational tool to track complexity variations of biosignals. FD analysis used in this work is based on Higuchi’s algorithm [15]. The algorithm is based on the measure of the mean length of the curve L(k) by using a segment of ksamples as a unit of measure.

Consider x(1), x(2),…, x(N) be the time series to be analysed. The algorithm constructs knew time series xmk, defined as

xmk:xm,xm+k,xm+2k,..……,xm+N−mkkform=1,2,….,kE14

where mand kare integers indicating the initial time and the interval time respectively. ameans the integer part of a.For each of the curves or time series xmkconstructed, the average length Lmkis computed as

Lmk=∑i=1N−mkxm+ik−xm+i−1kN−1N−mkkE15

where Nis the total length of the data sequence xand N−1N−mkkis a normalisation factor.

An average length is computed for all time series having the same delay (or scale) k,as the mean of the klengths Lmkfor m = 1,....,k. This procedure is repeated for each kranging from 1 to k_max, yielding a sum of average lengths L(k) for each kas indicated in Eq. (16).

Lk=∑m=1kLmkE16

The total average length for scale k, L(k) is proportional to k^−D, where Dis the FD by Higuchi’s method. In the curve of ln(L(k)) versus ln(1/k), the slope of the least squares linear best fit is the estimate of the FD. The k_maxvalues depends on the dimension D, on the signal’s length Nand on the specific class of fractal signals. The value of k_maxchosen for the analysis of EEG signals is 6. Figure 12 shows the plot of Fractal Dimension values calculated from the EEG signals of normal and depression patients acquired from the left part of the brain under eyes closed condition.

FD based on the algorithm followed is a quantifier evaluated directly in the time domain. Similar to the plot in Figure 12, the plots of FD for the rest of the recording protocols show that the values of FD are higher for normal controls indicating higher complexity in EEG signals of normal controls [16, 17]. It may be concluded that the value of FD increases with the increase in the degree of the cognitive activity. Lower values of FD indicate a low degree of cognitive activity for depression patients.

2.8.3 Largest Lyapunov Exponent (LLE)

Lyapunov exponents provide a qualitative and quantitative characterisation of dynamical behaviour. To discriminate between chaotic dynamics and periodic signals, Largest Lyapunov Exponent (λ) is often used. Lyapunov exponents are the average exponential rates of divergence or convergence of nearby orbits in phase space. It is a quantitative measure of the sensitive dependence on the initial conditions. The rate of separation between the nearby orbits in phase space is characterised by the Largest Lyapunov Exponent λ,mathematically written as:

Xτ−Yτ≈X0−Y0eλτE17

where X0and Y0are two initial conditions close together, and Xτand Yτare their respective time evolutions after τtime units. The approach of Wolf et al. [18] was used for constructing the algorithm for calculating the Largest Lyapunov Exponent (LLE) in this work.

Let x0→tdenote a reference trajectory passing through x0→0at time t = 0 and let x1→tdenote a trajectory passing through x1→0at time t = 0. The Largest Lyapunov Exponent λx0→is defined with respect to the reference orbit x0→by

λx0→=limt→∞limΔx0→→01tlogΔxt→Δx0→E18

where Δx0→is the Euclidean distance between the trajectories x0→tand x1→tat an initial time t = 0 and Δxt→is the Euclidean distance between the trajectories x0→tand x1→tat a later time t. In this definition x1→tcan be any trajectory that is initially infinitesimally close to x0→0at time t = 0. The correspondence between sensitivity to initial conditions and a positive Lyapunov exponent is obvious in equation number 17. An embedding dimension of 10 and a delay of 1 were used for calculating LLE. Figure 13 shows the values of LLE calculated for normal and depression patients from the left half of the brain under eyes closed condition.

Larger values of LLE observed for all the recording protocols for normal controls are indicative of higher brain activity [16, 17]. Therefore LLE can be effectively used for discriminating the EEG signals of normal controls and depression patients.

2.9 Surrogate data analysis

The method of using surrogate data in time-series analysis was introduced by Theiler et al. to validate that a given time-series is nonlinear. Nonlinear indices such as FD andLLE are computed for several surrogate data series. Their values are compared with that of the nonlinear index computed for the original data. The lack of any statistically significant difference is interpreted as the deviation from a linear process. Surrogate data is constructed by phase randomising the original time series and has the same linear features like mean, variance, histogram and power spectrum as the original data. This method of generating surrogate data is based on the amplitude adjusted Fourier transform method, which yields the same distribution of amplitudes but randomises the phases from the spectral aspect [19]. Tentative surrogate data are obtained by inverse Fourier transform.

To test for a statistical significance of difference in FDand LLEbetween the original and the surrogate data, 10 surrogate data series were generated to match each original signal. Let LLE(D) be the LLEof the original data, and let LLE(Si) be the LLEof the 10 surrogate series (i = 1, …, 10). The mean and standard deviation, SD of LLE(Si) (i = 1, …, 10) are estimated as LLE(S)and SD(LLE(Si)). The statistical significance measure σthen is computed as follows:

It follows a Student ttest distribution with 9 degrees of freedom (∼t_9[1-α/2]). For α = 0.05, the critical value of tis 2.26. Accordingly, when the σ> 2.26, the null hypothesis is rejected at the 5% probability level, and the original data are considered to contain nonlinear features.

Nonlinear indices (FD and LLE) are computed for several surrogate data series and their values are compared with the ones computed for the original series. The demonstration of significant difference in nonlinear indices between the original and surrogate data is supportive of the presence of nonlinearity in the original data. Tables 1 and 2 show the calculation of the statistical significance measure for both the normal controls and depression cases for the nonlinear indices, FD and LLE. The results prove that the original data contain nonlinear features since the statistical significance measure is greater than 2.26.