2nd Reading August 22, 2014 12:48 1440008

International Journal of Neural Systems, Vol. 24, No. 0 (2014) 1440008 (15 pages) c World Scientific Publishing Company  DOI: 10.1142/S0129065714400085

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

A PREDICTIVE MODELING APPROACH TO ANALYZE DATA IN EEG–fMRI EXPERIMENTS SAIDEH FERDOWSI Faculty of Electrical Engineering Shahrood University Shahrood 3619995161, Iran [email protected] SAEID SANEI Department of Computing Faculty of Engineering and Physical Sciences University of Surrey, Guildford, GU2 7XH, UK [email protected] VAHID ABOLGHASEMI Faculty of Electrical Engineering Shahrood University Shahrood 3619995161, Iran [email protected] Accepted 26 July 2014 Published Online 26 August 2014 In this paper, a novel technique based on blind source extraction (BSE) using linear prediction is proposed to extract rolandic beta rhythm from electroencephalogram (EEG) recorded in a simultaneous EEG–fMRI experiment. We call this method CLP–BSE standing for constrained-linear-prediction BSE. Extracting event-related oscillations is a crucial task due to nonphase-locked nature and inter-trial variability of this event. The main objective of this work is to extract rolandic beta rhythm to measure event-related synchronization (ERS) with acceptable signal-to-noise ratio (SNR). The extracted rhythm is utilized for constructing a regressor to analyze functional magnetic resonance imaging (fMRI). The proposed method is a semi-blind technique which uses a spatio-temporal constraint for beta rhythm extraction. This constraint is derived from recorded EEG signals based on the prior knowledge about the frequency and location of the source of interest. The main reason of employing linear prediction as an effective algorithm to extract the EEG rhythm is the ability of extracting sources which have specific temporal structure. Performance of the proposed method is evaluated using both synthetic and real EEG data. The obtained results show that the proposed technique is able to extract ERS effectively. The maximum percentage of ERS obtained by filtering is 152% while the obtained ERS by CLP–BSE is 214%. In another experiment, the extracted event-related oscillations in beta band are used to make the necessary regressor for fMRI analysis. The results of EEG–fMRI coregistration confirm that there are correlation between the extracted rolandic beta rhythm and simultaneously recorded fMRI. This conclude that, the results of EEG–fMRI combination support the reliability of CLP–BSE output. Keywords: Simultaneous EEG–fMRI; linear prediction; rolandic beta rhythm; blind source extraction.

1440008-1

2nd Reading August 22, 2014 12:48 1440008

S. Ferdowsi, S. Sanei & V. Abolghasemi

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

1.

Introduction

Electroencephalogram (EEG) and functional magnetic resonance imaging (fMRI) are noninvasive techniques which are increasingly used to study human brain functions. fMRI measures changes of blood oxygen level dependent (BOLD) in the brain and provides a high resolution 3D image showing the activate regions.1 EEG represents the neural activity inside the brain by measuring the changes of potentials over the scalp. In contrast to fMRI, EEG provides low spatial resolution and high temporal resolution information about brain functionality. Fusion of EEG and fMRI can therefore compensate for the existing shortcomings of each technique and provides a more comprehensive understanding of the neural activities. In spite of these advantages, simultaneous recording of EEG and fMRI poses some difficulties. EEG signals acquired in the magnetic field suffer mainly from gradient and ballistocardiogram (BCG) artifacts.2,3 Artifact removal is an essential step to obtain meaningful results in EEG–fMRI fusion. After removing the artifacts, a suitable technique is needed to extract meaningful features such as event-related potentials (ERPs) or event-related oscillations from EEG signal. In this research, we aim to develop an efficient technique to detect the changes of neural activity presented in the eventrelated oscillations. However, EEG signal recorded inside MRI scanner does not have quality as clean as signals recorded outside the scanner. Therefore, an effective technique is essential to improve the feature extraction procedure. Neural activity recorded by EEG oscillates in various frequency bands including delta (up to 4 Hz), theta (4–7 Hz), alpha (8–12 Hz), beta (13– 30 Hz), and gamma (> 30 Hz).4–6 The changes in synchronous activity of the neuronal population cause event-related changes in the amplitude of EEG oscillations.7 There are two strong event-related changes in brain oscillations: Event-related desynchronization (ERD) in alpha and beta bands and event-related synchronization (ERS) in beta band. Movement preparation suppresses EEG oscillation in both alpha and beta bands over the sensory motor area.8–10 This phenomenon is known as ERD and starts approximately 1 s before the movement. At the movement onset alpha amplitude decreases and

beta amplitude increases.4,11 This sharp power of beta band is known as ERS or post movement beta rebound (PMBR). Using the information obtained from EEG rhythms as a regressor for fMRI analysis gives the ability to localize the BOLD which is correlated with a specific neuronal rhythm. Most of the methods used for calculating ERD/ERS, measure the power changes of EEG signals in the corresponding frequency band obtained by averaging over a number of trials containing the movement. Since the event-related oscillation is a nonphase-locked response and varies from one trial to another, the extraction and analysis of ERD/ERS is a difficult task for single trial measurement. Studying variability in amplitude, latency or scalp distribution of ERD/ERS over consecutive trials reveals important information about cognitive and physiological states.12 Therefore, a method that can extract ERD/ERS in single trial recording with acceptable accuracy is essential for this purpose. Particularly, study on correlation between neural oscillations and BOLD requires time series in which ERD or ERS are demonstrated accurately. These time series are then used to make regressor for fMRI analysis techniques. Moreover, such a method allows us to have shorter experiment time which is helpful for subjects with impairment in motor or cognitive performance. Some techniques have been developed to analyze evoked brain potentials (time and phase locked response) and event-related oscillations (time and nonphase locked response) in single trial. Majority of approaches used for this purpose are based on blind source seperation (BSS). For example, in Ref. 13, independent component analysis (ICA) is used to identify and localize ERP in single trial multichannel EEG data. In another research, Cong et al.14 developed a non-negative multi-way array decomposition for multi-domain feature extraction for small ERPs. Lee et al.15 proposed an ICA-based spatio-temporal method to analyze single trial post-movement beta synchronization in magnetoencephalography (MEG) data. In another research, Ritter et al.16 used BSS to separate the rolandic alpha and beta rhythm from EEG signal recorded while performing a motor task. Their results show that the extracted sources i.e. rolandic alpha and beta rhythms are more correlated with motor task than raw EEG signals recorded at channels C3 and C4.

1440008-2

2nd Reading August 22, 2014 12:48 1440008

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

A Predictive Modeling Approach to Analyze Data

In this paper, we propose a blind source extraction (BSE) algorithm using linear predictor to extract beta rhythm from EEG signals in time domain. EEG dataset used in this work has been recorded in a joint EEG–fMRI experiment. In all the methods based on BSS, the recorded EEG signals are considered as weighted linear combinations of underlying cortical source signals. These weights are related to the distance between the cortical source regions or patches and electrodes, the orientation of the cortical patch relative to the electrodes, and the electrical properties of the intervening tissues (cortex, cerebral spinal fluid, skull, and skin). However, the inherent scaling ambiguity does not allow an exact estimations of these weights. To solve this problem several attempts have been made to incorporate the available signal characteristics and to develop semi-BSE techniques. In this work, the spatio-temporal information about the source of interest representing the EEG activity in beta band is used to develop a BSE algorithm to extract rolandic beta rhythm. This information is added as a constraint to a linear predictor. The constraint is called spatio-temporal since it is built up as a result of filtering the recoded EEG signal from a specific region of scalp. The main objective of the proposed method is to enhance the signal-to-noise ratio (SNR) for the extracted rhythm leading to better demonstrating of variations in the amplitude of this rhythm. Although, the proposed method is applied to the EEG signal in time domain, it also considers the information about the frequency content of the source of interest. This information can be achieved by bandpass filtering of signal selected as constraint around a specific frequency band. Computing PMBR requires a technique to obtain the power of extracted rolandic beta rhythm in different frequency bins. As a more efficient option for time-frequency analysis, compared to short-time Fourier transform (STFT), Adeli et al.17,18 used wavelet to analyze the EEG signals. In another research, Ahmadlou et al.19 proposed a spatiotemporal wavelet-chaos methodology to analyze EEG sub-bands for detecting abnormality in Alzheimer disease (AD). Here, we use complex Morelet wavelet to study the behavior of power of extracted rolandic beta rhythm. The computed wavelet energy of extracted rolandic beta rhythm is then used to make the regressors for general linear

model (GLM) in order to analyze the fMRI data. fMRI analysis is performed to assess the reliability of the output of proposed method. In the next section, the fundamental theory of BSS and BSE is briefly explained. In Sec. 3, the proposed method for extracting rolandic beta rhythms is described. The experimental results are presented in Sec. 4. In Sec. 5, the extracted rolandic beta rhythm is used to make a regressor for fMRI analysis. Finally, the paper is concluded in Sec. 6. 2.

Problem Formulation

What follows describes the BSS problem. Assume that x(t) = [x1 (t), x2 (t), . . . , xm (t)]T denotes m observed signals at time t, which are weighted sum of n individual sources s(t) = [s1 (t), s2 (t), . . . , sn (t)]T at time t (the superscript T represents vector transpose operator): xi (t) =

n 

aij sj (t) + vi (t)

i = 1, 2, . . . , m

(1)

j=1

or in matrix notation x(t) = As(t) + v(t),

(2)

where A denotes an m × n matrix called mixing matrix. In our application, columns of mixing matrix represent the projections from sources to the scalp electrodes. v(t) is an m × 1 noise vector. In BSS techniques, the source signals are retrieved by determining an n × m separation matrix W. In general, W can be represented by product of A−1 and scaling and permutation matrices. After estimating the separation matrix, the source signals are calculated by: yj (t) =

m 

wji xi (t)

j = 1, 2, . . . , n

(3)

i=1

or in matrix notation y(t) = Wx(t),

(4)

where y(t) = [y1 (t), y2 (t), . . . , yn (t)]T is the n × 1 estimated source signals at time t. In contrast to BSS algorithms which recover all the sources simultaneously, BSE algorithms recover only one or a selected number of sources each time. The problem of estimating sources of interest in BSE is equivalent to identification of the corresponding columns of mixing matrix, aj , or the rows of its

1440008-3

2nd Reading August 22, 2014 12:48 1440008

S. Ferdowsi, S. Sanei & V. Abolghasemi

pseudo inverse, wj . Extracting one source at each time, enables us to impose constraints while the algorithm is estimating the source of interest. A BSE algorithm may exhibit weaker convergence than simultaneous blind separation techniques. The main reason is the existence of some illconditioned problems due to accumulation of error during deflation procedures. In order to eliminate this problem a pre-withening is applied to the observation signals x in the form of

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

˜ (t) = Qx(t), x

(5)

where Q ∈ Rm×m is a decorrelation matrix ensuring ˜ T } = Im , x·x that the correlation matrix, Rx˜ = E{˜ is an identity matrix. It has been shown that a BSE approach using linear predictor is capable of extracting temporally correlated sources.20 Since this is in the scope of this paper, linear prediction is described in detail in the following section. 3.

Source Extraction Using Linear Prediction

In source extraction using linear prediction, the predictability of source signals is used as the main criterion for estimating the rows of separation matrix, denoted by wj . Based on this assumption, yj (t) can be modeled/predicted using a weighted linear combination of previous samples.a y(t) = e(t) +

K 

T

¯ (t), bp y(t − p) = e(t) + b y

(6)

p=1

where e(t) represents the prediction error, K is prediction order, b = [b1 , b2 , . . . , bK ]T are weights of ¯ (t) = [y(t − 1), y(t − 2), . . . , y(t − prediction and y K)]T . Equation (6) is written as follows by substituting y(t) with wT x(t): e(t) = wT x(t) −

K 

bp wT x(t − p),

where w = [w1 , w2 , . . . , wm ] . It is desired to estimate the optimal value of w and b leading to a successfully extracted y. For this purpose, the prediction error should be minimized: 1 (8) J (w, b) = E{e2 (t)}, 2 a

3.1. Proposed method to extract rolandic beta rhythm In this section, details of the proposed method for extracting rolandic beta rhythm is presented. This method extends the model introduced in linear prediction Eq. (8) to exploit the prior information about the location, temporal structure, and bandwidth of the source of interest. This information are imposed using a reference signal obtained by bandpass filtering of recorded EEG signal by electrodes placed in the location expected to observe source of interest. What follows presents the proposed constrained minimization problem: min J (w, b), w,b

s.t. min(ξ − E{(y(t)r(t))2 })2 ,

(9)

where r = [r(1), r(2), . . . , r(T )]T is the reference signal defined based on the location and frequency content of the source of interest and ξ is a suitably chosen positive constant. The above constrained minimization problem changes to an unconstrained function as follows: JCLP (w, b) =

1 E{e2 (t)} 2 + λ(ξ − E{(y(t)r(t))2 })2 ,

(7)

p=1 T

where E{·} denotes the expected value. Minimizing the above cost function leads to extracting sources with specific temporal structure. However, the performance could be significantly improved by imposing a priori about the temporal structure of sources of interest. This can play a crucial role in applications concerning artifact removal from biomedical signals, especially EEG signals. The obtained results presented in Sec. 4 demonstrate the advantage of using the priori knowledge during source extraction procedure.

(10)

where λ is the regularization parameter which stabilizes the tradeoff-ratio between J (w, b) and the constraint. The standard stochastic gradient descent algorithm21 is used to estimate w and b which minimizes the proposed cost function. The partial derivatives of JCLP (w, b) with respect to w and b are calculated

For notational simplicity, we omit index j from yj (t) and wj in the sequel. 1440008-4

2nd Reading August 22, 2014 12:48 1440008

A Predictive Modeling Approach to Analyze Data

as follows:

Algorithm 1: CLP-BSE.

∂JCLP (w, b) = wT − 2wT ∂w + wT

K 

Input: X: observation matrix, K: short-term prediction orders, λ0 : initial value for regularization parameter, γ: positive scaler within (0, 1] and n: number of sources. Output: W: separation matrix. begin for j = 1 to n do random initialization of vectors wj , b. repeat ∂J (w,b) • wj ← wj − ηw ( CLP ) ∂w • wj ← wj /wj 2 • ∀ k = 1, . . . , K : ∂J (w,b) bk ← bk − ηb ( CLP ) ∂bk

bp Rx˜ (p)

p=1

K  K 

bp bq Rx˜ (p − q)

p=1 q=1

− 4λ(ξ − wT Rx˜r (0)w) × (wT Rx˜r (0)),

(11)

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

∀ k = 1, . . . , K: ∂JCLP (w, b) = −E{e(t)y(t − k)} ∂bk = −wT Rx˜ (k)w +

K 

1−γ • λ ← λ0 1−γ itr • itr ← itr + 1

T

bp w Rx˜ (p − k)w,

until a stopping criterion is met ; end end

(12)

p=1

where bk is kth element of b. The learning rules for w and b are derived by substituting the derivatives obtained in (11) and (12) into the stochastic gradient descent approach.  K T T  bp Rx˜ (p) w(u+1) = w(u) − ηw w(u) − 2w(u) p=1

+ w(u)

T

K  K 

bp bq Rx˜ (p − q)

p=1 q=1

− 4λ(ξ − w

(u)T



Rx˜r (0)w)(w

(u)T

Rx˜r (0)) , (13)

∀ k = 1, . . . , K: (u+1)

bk

1−γ , (15) 1 − γ itr where λ0 is initial value of λ, γ is a positive scalar selected within (0 1] and itr is the iteration counter. The above updates are carried out for all the sources from 1 to n and the procedure of updating the values of w and b is repeated to minimize (10). It is important to note that the normalization of w is necessary after each iteration to preserve the column norm of separation matrix. λ = λ0

4.  T

(u)

= bk − ηb − w(u) Rx˜ (k)w(u)

+

variation:

K 

 T

(u) b(u) Rx˜ (p − k)w(u) , p w

(14)

p=1

where ηw and ηb are step sizes and w(u+1) and (u+1) are the updated values of parameters of the bk proposed cost function. We used fixed step sizes for this method. A pseudo-code of CLP–BSE is given in Algorithm 1. The regularization parameter λ, is made variable starting from a high value and decreasing toward zero. We used the following function for this

Experimental Results

Two sets of data including synthetic and real EEG data are used to evaluate the performance of the proposed methods. The results are represented in this section. Then, the extracted source of interest from real EEG dataset which demonstrates the ERS in beta band is used to build a regressor for fMRI analysis. For this purpose, we use Morlet wavelet transform to detect the instantaneous interaction between the power of EEG signal in the beta band and motor task. 4.1. Synthetic data Ten sources including EEG rhythms and noise are used to generate a set of synthetic EEG signals

1440008-5

2nd Reading August 22, 2014 12:48 1440008

S. Ferdowsi, S. Sanei & V. Abolghasemi

0.2 s1

0 −0.2 0 0.2

s2

Amplitude

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

1

2

3

4

5

theta

upper beta 1

2

3

4

5

1

2

3

4

5

1

2

3

4

5

1

2

3

4

5

1

2

3

4

5

1

2

3

4

5

1

2

3

4

5

2

3

4

5

lower beta

noise delta

noise

gamma

alpha

Fig. 2.

Simulated source locations inside the brain.

0

−0.2 0 0.2 s8 0 −0.2 0 0.2 s9 0 −0.2 0 0.2 s10 0 −0.2 0

1

Time (second)

Fig. 1.

Simulated EEG sources.

(Fig. 1). The EEG rhythms are derived from real EEG data. These rhythms belong to five different frequency bands: delta (s1), theta (s2), alpha (s3), beta (s4 and s5), and gamma (s6). s4 is related to lower beta band (13–17 Hz), s5 represents the upper beta band (17–26 Hz), and s7–s10 are random uniform noise. The sampling rate of synthetic EEG signals is 128 Hz. In order to generate the synthetic EEG signals similar to the real scenario, we followed the model of brain described in Refs. 22–24. In this model the EEG signals denoted by X in (16) are represented by an m × T matrix, where m is the number of electrodes and T is the number of time samples. X = HMS + N =

m 

Hq mq sq + N.

(16)

q=1

In this expression, H is an m × 3n matrix describing the forward mixing model of the n sources to the m electrodes, M is a 3n × n matrix describing the orientation of the n sources and S is an n × T matrix containing the time-course of sources. In this work, BrainStorm softwareb is used to create the forward b

noise

0 −0.2 0 0.2

s7

5

0 −0.2 0 0.2

s6

4

0 −0.2 0.20

s5

3

0 −0.1 0 0.2

s4

2

0 −0.2 0 0.1

s3

noise 1

model. We used a three layer spherical head model with conductivities of 0.33, 0.0042, and 0.33 µS/cm, for scalp, skull, and brain, respectively. The generated synthetic sources are placed in different locations inside the brain. Source locations are shown in Fig. 2. Then, EEG signals in 19 channels are obtained using (16). In the last step, Gaussian noise is added to all the channels to obtain synthetic data with different SNRs. The reference signal r is obtained based on the knowledge about the location of the source of interest. In this work, we are interested to extract EEG oscillations in upper-beta band in motor cortex. So, a (17–25 Hz) bandpass filtering in channel C3 is undertaken to obtain the reference signal. Prior to applying the proposed method to synthetic EEG signal, we need to find the proper value of prediction order (p). Therefor in the first experiment, the algorithm is executed for different values of prediction orders started from 2 to 25 to find the optimal value of prediction order. It should be noted that, in this experiment λ is set to zero. The following performance index is used to evaluate the performance of the algorithm for different prediction order.21 In all source extraction algorithms, it is desired to find a vector wj such that the value of each source at time t is estimated using yj (t) = wjT x(t). In the ideal case, gj = AT wj = ±cj ej is a vector with only one nonzero element in the jth place, corresponding to the jth extracted source. ej is a vector with only one unitary element, and cj is a positive scalar. Vector gj is used to quantify the quality of source extraction algorithm as defined in the following

http://neuroimage.usc.edu/brainstorm/. 1440008-6

2nd Reading August 22, 2014 12:48 1440008

A Predictive Modeling Approach to Analyze Data

recorded data channels (in kbits/sec).25 The amount of mutual information removed from the results of separation is calculated as follows:

−8 −10

PI 3 = MIR = I(x) − I(y) = [h(x1 ) + · · · + h(xm )]

−14

1

PI (dB)

−12

− [h(y1 ) + · · · + h(yn )] + log |W|,

−16

(19)

−18 −20 −22

0

5

10

15

20

25

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

p

Fig. 3.

Numerical values of PI 1 versus p.

equation: PI 1 =

   n n 2  gij 1 1 , 10 log10   2 , . . . , g2 } − 1 n i=1 n j=1 max{gi1 in (17)

where n is the number of extracted sources. Smaller value for this index indicates better quality in the results of source extraction. Figure 3 shows the obtained values for PI 1 when prediction order varies between 2 and 25. It is seen that this index reaches to its minimum value for 10 ≤ p ≤ 17. Based on the obtained results for this experiment the prediction order is set to 12 in the later simulations. After finding the best value for prediction order, in the second experiment CLP–BSE is applied to synthetic dataset in order to extract source of interest. In this step, the performance of the proposed algorithm is compared with other BSS methods. Two performance indices (PI 2 and PI 3 ) are used to show the capability of the proposed method to extract the source of interest. PI 2 is a performance index which measures the correlation between the extracted source as the beta rhythm and the actual one. This index is calculated as follows.     y(t), s(t)   − 1, (18) PI 2 = 10 log10   y(t), s(t)y(t), s(t) 

where I(x) and I(y) are the mutual information of recorded signals and recovered sources respectively. h(xi ) and h(yj ) represent respectively the marginal entropy of vectors x and y. Tables 1 and 2 present the values of PI 2 and PI 3 for LP–BSE, CLP–BSE, Infomax, and SOBI when different amounts of noise are added to the mixture. LP–BSE indicates the source extraction based on linear prediction when no constraint is incorporated. Infomax and SOBI are well known ICA algorithms proposed by Bell and Sejnowsky26 and Belouchrani et al.27 respectively. As can be seen from Table 1, the proposed method outperforms other BSS techniques specially when more noise is added to the mixture. In addition to imposing the constraint, the assumption of having colored sources by auto-regressive model leads to improving the results. This assumption makes the linear predictor-based techniques suitable for extracting periodic and semi-periodic sources Table 1. Obtained values for PI 2 (in dB) for different methods. Method SNR

LP–BSE

CLP–BSE

Infomax

Sobi

no-noise 20 15 10 5

−19.46 −15.49 −12.93 −9.82 −7.59

−26.96 −22.43 −16.53 −12.63 −11.34

−16.81 −14.08 −10.84 −8.29 −7.03

−22.29 −18.18 −13.15 −10.17 −8.65

where y(t) is the source of interest extracted by the algorithm and s(t) is actual source. Smaller value of this index show better performance of the algorithm. The third performance index, used for comparison, measures the mutual information reduction (MIR) between the recovered sources relative to the 1440008-7

Table 2. Obtained values for MIR (in kbits/s) for different methods. Method SNR

LP–BSE

CLP–BSE

Infomax

Sobi

no-noise 20 15 10 5

37.03 32.97 29.77 29.81 25.61

37.18 33.55 31.06 29.75 27.03

37.14 33.08 30.59 28.02 25.74

37.13 33.24 30.92 29.31 26.49

2nd Reading August 22, 2014 12:48 1440008

10

Original source LP−BSE CLP−BSE Infomax Sobi

Amplitude

such as EEG rhythms. It is also seen that SOBI which is based on second order statistics and uses similar assumption about the sources, shows highest performance after CLP–BSE. The results presented in Table 2 show the obtained values for MIR by different algorithms. It is seen that all the methods provide nearly the same MIR for noiseless data which is about 37 kbits/s. It is seen that Informax, SOBI, and LP–BSE have poorer performance than CLP–BSE in removing mutual information when the level of noise increases. It seems that incorporating constraint improves this index, though not considerable. The above results obtained, when 0.5 and 0.9 are selected for parameters λ0 and γ respectively. Figure 4 represents the trend of changing values of λ during the source extraction procedure. As seen, the effect of constraint is higher at the beginning but it decreases as the iterations evolved. The extracted source of interest recovered by different methods, when the input SNR is 5 dB, are presented in Fig. 5(a). As it is seen from the figure, the extracted beta rhythm using CLP–BSE well approximates the original source. The estimated weights corresponding to the extracted beta rhythm are also shown in Figs. 5(b)–5(e). As it is seen from Figs. 5(b)–5(e), the topographic maps after projection show that the maximum weight refers to the left side of motor area where the original source in upper beta rhythm has been placed. Comparing the topographic maps reveals that the estimated weights by infomax (Fig. 5(d)) do not localize the correct origin of upper beta activity. Moreover, the Talariach coordinates of the recovered sources and dipolarity are estimated in this experiment. Dipolarity of the source extraction is defined

0

1

2

Time (second)

3

4

5

(a)

(b) LP–BSE

(c) CLP–BSE

(d) Infomax

0

(e) SOBI 10

Fig. 5. (a) Extracted source of interest, and (b)–(e) topographic projection of the estimated mixing weights for source of interest, when the input SNR has been 5 dB.

−1

λ

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

S. Ferdowsi, S. Sanei & V. Abolghasemi

10

10

−2

−3

0

50

100

150

200

250

300

350

400

450

500

Iteration

Fig. 4.

Trend of updating λ versus number of iterations.

as the number of estimated components which their scalp maps can be fit to the scalp projection of a single equivalent dipole with an error less than a specified error threshold.25 Table 3 presents the

1440008-8

2nd Reading August 22, 2014 12:48 1440008

A Predictive Modeling Approach to Analyze Data

Table 3. Talairach coordinates and dipolarity of recovered source of interest by different methods when the input SNR has been 5 dB. Method

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

LP–BSE CLP–BSE Infomax Sobi

Talairach coordinate (x, y, z)

ND%

(−50, 14, 48) (−50, 2, 47) (−18, 68, 18) (−57, −8, 55)

33.81 38.84 32.45 34.21

Talariach coordinates and the percentage of dipolarity for the extracted source of interest when noise level is 5 dB. It should be noted that the exact coordinate of upper beta rhythm is (−41.28, −2.67, 47.24). Comparing the results reveal that the proposed method provides more accurate results than the other techniques. 4.1.1. Simultaneous EEG–fMRI recording The real EEG data used in this study were recorded from five healthy men. All the subjects were righthanded and their ages ranged from 18 to 50 years. The study was approved by the local ethics committee of King’s College London. The EEG was acquired using the Neuroscan Maglink RT system (impedance was kept within 10–20 Kohms), providing 64-channel comprising 62 scalp electrodes, one ECG electrode, and one EOG electrode. The sampling rate of raw EEG data was set at 10 KHz. During the fMRI acquisition, 300 volumes including 38 slices (3.2969 × 3.2969 × 3.3 mm resolution, TR = 2000 ms, TE = 25 ms) were acquired. In order to decrease the computational time, the number of channels decreased to 32 including O2, O1, OZ, PZ, P4, CP4, P8, C4, TP8, T8, P7, P3, CP3, CZ, FC4, FT8, TP7, C3, FZ, F4, F8, T7, FT7, FC3, F3, FP2, F7, FP1, PO3, PO4, PO8, PO7. All the preprocessing steps used for this dataset were: (1) removing gradient or imaging artifact using the method proposed by Niazy et al.;28 (2) downsampling the data to 250 Hz; (3) bandpass filtering using a Butterworth filter between 0.5 Hz and 45 Hz; (4) in the last step, reference free recordings were obtained using the common average referencing; (5) removing BCG artifact using the method developed by Ferdowsi et al.3 The experimental design consists of three blocks: Cued movement, free movement, and visual control. During the cued movement, one of the arms of a cross

would become illuminated (green) and participants were required to move the joystick in the illuminated direction, during the free movement two arms would become illuminated (red), and participants were free to move toward one of the nonilluminated arms and during the visual control block one of the arms of a cross changed color (red) and participants should not move the joystick. In all three blocks, there were nine trials per block and each trial including light turning on and then off lasted for 2.4 s. The proposed method is applied to the segments of real EEG signal obtained after removing BCG artifact using method proposed in Ref. 3. Although output of BCG removal algorithm is considered as a clean signal, it still suffers from other contaminating biological and measurement induced interferences such as residual of gradient artifact. The main aim of CLP–BSE, is to extract rolandic beta rhythm from such a dataset. In this section, the results of applying CLP–BSE to real EEG signals are presented. Figure 6 shows the results of factorization for one sample segment of data. The third extracted source represents the brain activity in beta band. The extracted topographic map also shows the region of activity in the left area of primary motor cortex. In order to obtain the time-frequency transform of the signal we use wavelet transform. Complex Morlet wavelet is a useful technique to provide time-varying representation of the energy of the EEG signal. The complex Morlet wavelet w(t, f0 ) is formulated as follows29 :

2 −t (20) w(t, f0 ) = A exp exp(2iπf0 t) 2σt2 √ with σf = 1/2πσt and A = (σt π)1/2 . The tradeoff ratio (f0 /σf ) is equal to 7 to create a wavelet family. We created a family of Morlet wavelet at 1 Hz frequency intervals in the range of 13 to 30 Hz. The time-varying energy E(t, f0 ) of a signal at a specific frequency band can be calculated by convolving the complex wavelet with signal (Eq. (21)). E(t, f0 ) = |w(t, f0 ) ∗ y(t)|2 .

(21)

Figure 7(a) shows the estimated power timecourse of y3 . The calculated time-frequency transform for rolandic beta rhythm obtained by filtering the EEG signal from channel C3 is also shown in Fig. 7(b). Pink dashed lines in the figure show the movement triggers. In is seen that, PMBR in Fig. 7(a) is more visible than in Fig. 7(b).

1440008-9

2nd Reading August 22, 2014 12:48 1440008

S. Ferdowsi, S. Sanei & V. Abolghasemi

y1 y2 y3 y4 y5

Amplitude

y6 y7 y8 y9 y10 y11 y12 y13 y14 y15 y16 y17 y18

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

y19 y20 0

1

2

3

4 Time(second)

5

6

7

8

(a)

(b)

(c)

Fig. 6. Extraction results using CLP–BSE; (a) The estimated source signals and (b)–(c) Right and left view of topographic maps of the estimated sources. y3 , as the extracted rolandic beta rhythm, presents brain activity in motor cortex area.

(a)

Fig. 7.

(b)

Power spectrum of the rolandic beta rhythm; (a) estimated by CLP–BSE and (b) extracted by filtering.

1440008-10

2nd Reading August 22, 2014 12:48 1440008

A Predictive Modeling Approach to Analyze Data

Table 4.

Percentage of calculated ERS for different techniques (%).

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

12–16 Hz

17–22 Hz

23–30 Hz

Subject

Filtering

CLP–BSE

Filtering

CLP–BSE

Filtering

CLP–BSE

1 2 3 4 5

73 67 62 53 47

96 89 87 85 69

152 183 125 96 89

214 236 198 239 172

25 46 45 34 27

38 59 67 47 41

In order to evaluate the results of the proposed method, the ERS values are calculated for the results of the proposed method and convectional filtering technique. A band-pass butterworth filter with suitable cut-off frequencies is used to filter the EEG signal and then wavelet transform is used to calculate the power time-course. The recorded EEG by channel C3 is used as input of the filter. The following equation is used to calculate the ERS. ERS =

(Pactivation − Prest ) × 100, Prest

(22)

where Pactivation is the signal power when beta rebound occurs and Prest is the signal power during the specific reference interval. ERSs are calculated for three frequency bands covering the beta rhythm. Table 4 presents the averaged values for ERS of each subject. The obtained results for ERS show that the spatio-temporal constraint increase ERS values. Higher values for ERS resulted from improving the SNR of the signal representing the brain activity in beta band. The frequency bin which has the maximum power in each band is used to calculate the ERS. The highest values for ERS are obtained for frequencies between 17 Hz and 22 Hz. As described above, the power time-course of the EEG signal in beta band is obtained using wavelet transform. The calculated power time-course represents the instantaneous interaction between the EEG activity in beta band and motor task. Power time-course is then convolved with the Hemodynamic response function (HRF) to make the regressor for fMRI analysis. This regressor is used to predict the BOLD response in fMRI data which are collected simultaneously with EEG. A schematic illustration of data analysis is given in Fig. 8. The next section presents the results of fMRI analysis.

Fig. 8. Schematic representation of different steps of EEG–fMRI analysis.

1440008-11

2nd Reading August 22, 2014 12:48 1440008

S. Ferdowsi, S. Sanei & V. Abolghasemi

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

5.

fMRI Data Analysis

Preprocessing of functional MRI including motion correction in three dimension, normalization, and smoothing were carried out using SPM toolbox. Removing nonneural component, the so-called physiological noise (PN) is also considered as a further step to obtain less noisy results. The main cause of PN is variation in the subject’s heart rate and respiratory rate and arterial carbon dioxide.30 In this work, the PN is filtered using GLM.31–33 For this purpose, first a template of noise is generated by averaging the time courses of the voxels of brain ventricle. Then, this template is used by GLM as a regressor to filter nonneural noise from fMRI. After preprocessing, areas of significant BOLD contrast are determined using GLM. A box car design matrix is created using the stimulus onset timings. Motion correction parameters obtained during the preprocessing also is included to account for subject movements. A fixed effect group analysis is performed. Figure 9 shows the results of group analysis. Highlighted regions present the activated areas due to performing the motor task. Group activation

maps are thresholded at p < 0.05 (Bonferronicorrected) and are overlayed on the structural scan for a subject in Talairach space. The contrast images are calculated for two types of movements, Cued and Free versus Control. The colored voxels highlight the areas where the regressors (Cued and Free movements) accounted for significant variance in the MRI images. It is seen that the significant region is located in the primary motor cortex. This confirms the results obtained for mixing column corresponding to source of interest in previous section. Since the active area obtained for both types of movement are overlapped, the contrast images are calculated for both types of movements (Cued and Free) versus Control. It should be noted that removing the PN leads to increasing the T -value for active regions which means having less variance of the voxel signals. In this experiment, the T -value increased 2.08 unit after removing the PN noise. In another experiment, the correlation between the extracted power time-course and fMRI time series is identified using a separate model for each subject. Figure 10(a)–10(e) shows the correlation

(a) Cued movement

(b) Free movement

Fig. 9. Regions of motor activation over all subjects. In these figures SAG, COR, and TRA refer respectively to sagittal, coronal, and transversal views.

(a) Subject 1

(b) Subject 2

Fig. 10. Combining EEG and fMRI to investigate the correlation of BOLD and PMBR for (a)–(e) individual subjects and (f) over all subjects. In these figures SAG, COR, and TRA refer respectively to sagittal, coronal, and transversal views. 1440008-12

2nd Reading August 22, 2014 12:48 1440008

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

A Predictive Modeling Approach to Analyze Data

(c) Subject 3

(d) Subject 4

(e) Subject 5

(f) Group analysis

Fig. 10.

between the beta rhythm and BOLD signal for individual subjects. In addition to EEG–fMRI combination for individual subject, group analysis also performed to detect the regions related to PMBR over all five subjects. The result obtained for group analysis is presented in Fig. 10(f). Removing the PN in this experiment improves the obtained T -values for 0.72 in average. The active maps obtained based on the energy of the extracted signal in beta band are very similar to the maps obtained after fMRI analysis using box car design matrix. In other words, the positive correlation between BOLD and wavelet energy confirms the findings by fMRI block analysis. In order to make a better comparison the number of detected active voxels in the region of interest is computed. The number of voxels in active area in motor cortex are 73, 48, 78, 89, and 57 for subjects 1, 2, 3, 4, and 5, respectively. All five subjects have demonstrated the region of activity of positive correlation between the PMBR and MRI images in the sensory motor cortex. 6. Discussion and Conclusion Combining information obtained from EEG and fMRI leads to a comprehensive understanding of brain functionality. Most common techniques used

(Continued )

for EEG–fMRI integration, utilize the extracted information from EEG to make a regressor for fMRI analysis. Two main objectives are investigated in this research; First, developing an effective technique to detect post-movement beta rebound or ERS in beta band in EEG signals obtained from a joint EEG– fMRI experiment. Second, investigating the correlation between detected PMBR and ERS with fMRI. Since the EEG signals simultaneously recorded with fMRI are more noisy than the signals recorded outside the MRI scanner, an effective technique for the noise removal is essential. We addressed this problem by proposing a new method called CLP– BSE. The proposed technique is a BSE algorithm that uses linear prediction to extract the sources. CLP–BSE proposed a semi-blind technique based linear prediction which uses the a priori information about the location, bandwidth, and temporal structure of rolandic beta rhythm through the source extraction procedure. Linear predictor is a suitable technique to model the colored or periodic sources such as EEG brain rhythms. Moreover, imposing the proposed constraints leads to incorporating the additional information and improving the quality of the extracted source of interest. Although, imposing constraint can improve the results, using reference signal as constraint which does not have

1440008-13

2nd Reading August 22, 2014 12:48 1440008

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

S. Ferdowsi, S. Sanei & V. Abolghasemi

enough quality may deteriorate the performance of the algorithm. Therefore, the quality of reference signal which carries the prior information should be investigated to achieve an acceptable convergence. The efficiency of CLP–BSE is examined using both synthetic and real data. In order to simulate the condition of MRI scanner for synthetic data, different noise levels are added to the simulated EEG signals. Comparing the results of proposed method and other BSS algorithms shows the ability of proposed method to extract the beta rhythm with acceptable SNR. The superiority of CLP–BSE to other BSS algorithms is more visible when the noise level is high. Moreover, the Talariach coordinates of the extracted beta rhythm by CLP–BSE and the selected benchmarks are also estimated. Since for the synthetic data, the exact location of each source in the brain is known, investigating the performance of all the algorithms in terms of source localization is possible. It is seen that, the proposed method returns a closer coordinate to the exact location of the synthetic beta rhythm. The main reason of this success is because of imposing constraint containing the prior knowledge about location of the beta rhythm. The performance of the proposed method is also investigated using real EEG signal recorded simultaneously with fMRI. After applying the proposed method to real EEG signals, the event-related synchronization is calculated using extracted rolandic beta rhythm. The obtained values of ERS are then compared with the results of conventional filtering technique. It is seen that the proposed method improves the visibility of changes in the power of beta rhythm due to performing a task reflected in the EEG signals. In addition, the EEG–fMRI co-registration to detect the active areas correlated with the extracted rolandic beta rhythm is performed. EEG–fMRI analysis leads to detecting the active areas in motor cortex which are positively correlated with the post-movement beta rebound. The concept of using constrained BSE techniques for joint EEG–fMRI analysis is an open field of research. Developing constrained techniques which use more a priori knowledge about the source of interest can improve the results of extracting features from EEG and fusing EEG–fMRI consequently. As one of our future work, this concept can be used to detect the changes in other EEG

brain oscillations. Moreover, investigating correlation between the BOLD component in fMRI and ERD in alpha band detected by the proposed method can be addressed in future. References 1. D. Rangaprakash, X. Hu and G. Deshpande, Phase synchronization in brain networks derived from correlation between probabilities of recurrences in functional MRI data, Int. J. Neural Syst. 23(2) (2013) 1350003. 2. P. J. Allen, O. Josephs and R. Turner, A Method for removing imaging artifact from continuous EEG recorded during functional MRI, NeuroImage 12(2) (2000) 230–239. 3. S. Ferdowsi, S. Sanei, V. Abolghasemi, J. Nottage and O. O’Daly, Removing ballistocardiogram artifact from EEG using short- and long-term linear predictor, IEEE Trans. Biomed. Eng. 60(7) (2013) 1900–1911. 4. S. Sanei and J. A. Chambers, EEG Signal Processing (Wiley, Chichester, UK, 2007). 5. S. Sanei, Adaptive Processing of Brain Signals (Wiley, Chichester, UK, 2013). 6. S. Ghosh-Dastidar, H. Adeli and N. Dadmehr, Mixed-band wavelet-chaos-neural network methodology for epilepsy and epileptic seizure detection, IEEE Trans. Biomed. Eng. 54(9) (2007) 1545–1551. 7. G. Pfurtscheller and F. H. Lopes da Silva, Eventrelated EEG/MEG synchronization and desynchronization: Basic principles, Clin. Neurophysiol. 110(11) (1999) 1842–1857. 8. H. Jasper and W. Penfield, Electrocorticograms in man: Effect of voluntary movement upon the electrical activity of the precentral gyrus, Europ. Arch. Psych. Clin. Neuro. 183 (1949) 163–174. 9. R. de Jong, T. E. Gladwin and B. M. Hart, Movement-related EEG indices of preparation in task switching and motor control, Brain Res. 1105(1) (2006) 73–82. 10. D. J. Krusienski, G. Schalk, D. J. McFarland and J. R. Wolpaw, A µ-rhythm matched filter for continuous control of a brain-computer interface, IEEE Trans. Biomed. Eng. 54(2) (2007) 273–280. 11. G. Pfurtscheller, A. Stancak and G. Edlinger, On the existence of different types of central beta rhythms below 30 Hz, Electroencephalograp. Clin. Neurophysiol. 102(4) (1997) 316–325. 12. T. Jung, S. Makeig, M. Westerfield, J. Townsend, E. Courchesne and T. J. Sejnowski, Analysis and visualization of single-trial event-related potentials, Hum. Mapp. 14(3) (2001) 166–185. 13. J. Tzyy-Ping, M. Scott, W. Marissa, T. Jeanne, C. Eric and J. S. Terrence, Analysis and visualization of single-trial event-related potentials, Hum. Brain Mapp. 14(3) (2001) 166–185.

1440008-14

2nd Reading August 22, 2014 12:48 1440008

Int. J. Neur. Syst. Downloaded from www.worldscientific.com by CHINESE UNIVERSITY OF HONG KONG on 12/27/14. For personal use only.

A Predictive Modeling Approach to Analyze Data

14. F. Cong, A. H. Phan, P. Astikainen, Q. Zhao, Q. Wu, J. K. Hietanen, T. Ristaniemi, and A. Cichocki, Multi-domain feature extraction for small eventrelated potentials through nonnegative multi-way array decomposition from low dense array EEG, Int. J. Neural Syst. 23(2) (2013) 1350006. 15. P. Lee, Y. Wu, L. Chen, Y. Chen, C. Cheng, T. Yeh, L. Ho, M. Chang and J. Hsieh, ICA-based spatiotemporal approach for single-trial analysis of postmovement MEG beta synchronization? NeuroImage 20(4) (2003) 2010–2030. 16. P. Ritter, M. Moosmann and A. Villringer, Rolandic alpha and beta EEG rhythms strengths are inversely related to fMRI-BOLD signal in primary somatosensory and motor cortex, Hum. Brain Mapp. 30(11) (2009) 1168–1187. 17. H. Adeli and S. Ghosh-Dastidar, Automated EEGBased Diagnosis of Neurological Disorders: Inventing the Future of Neurology (Taylor & Francis, Boca Raton, FL, USA, 2010). 18. H. Adeli, S. Ghosh-Dastidar and N. Dadmehr, A spatio-temporal wavelet-chaos methodology for EEG-based diagnosis of alzheimer’s disease, Neurosci. Lett. 444(2) (2008) 190–194. 19. M. Ahmadlou, H. Adeli and A. Adeli, Fractality and a wavelet-chaos-methodology for eeg-based diagnosis of alzheimer disease, Alzheimer Dis. Assoc. Disord. 25(1) (2011) 85–92. 20. A. Cichocki and R. Thawonmas, On-line algorithm for blind signal extraction of arbitrarily distributed, but temporally correlated sources using second order statistics, Neural Process. Lett. 12(1) (2000) 91–98. 21. A. Cichocki and S. I. Amari, Adaptive Blind Signal and Image Processing: Learning Algorithms and Applications (John Wiley & Sons, Inc., New York, NY, USA, 2002). 22. P. Schimpf, J. Haueisen, C. Ramon and H. Nowak, Realistic computer modelling of electric and magnetic fields of human head and torso, Parallel Comput. 24(9–10) (1998) 1433–1460.

23. J. Sarvas, Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem, Phys. Med. Biol. 32(1) (1987) 11–22. 24. L. Spyrou and S. Sanei, Source localization of eventrelated potentials incorporating spatial notch filters, IEEE Trans. Biomed. Eng. 55(9) (2008) 2232–2239. 25. J. Onton, R. Oostenveld, S. Makeig, A. Delorme and J. Palmer, Independent EEG sources are dipolar, PLoS ONE 7(2) (2012) e30135. 26. A. J. Bell and T. J. Sejnowski, An informationmaximization approach to blind separation and blind deconvolution, Neural Comp. 7(6) (1995) 1129–1159. 27. A. Belouchrani, K. Abed-Meraim, J. F. Cardoso and E. Moulines, A blind source separation technique using second-order statistics, IEEE Trans. Signal Proces. 45(2) (2002) 434–444. 28. R. K. Niazy, C. F. Beckmann, G. D. Iannetti, J. M. Brady and S. M. Smith, Removal of fMRI environment artifacts from EEG data using optimal basis sets, Neuroimage 28(3) (2005) 720–737. 29. A. Teolis, Computational Signal Processing with Wavelets (Applied and Numerical Harmonic Analysis) (Birkhauser, Boston, USA, 1998). 30. R. M. Birn, K. Murphy, D. A. Handwerker and P. A. Bandettini, fMRI in the presence of taskcorrelated breathing variations, Neuroimage 47(3) (2009) 1092–1104. 31. K. J. Friston, A. P. Holmes, K. J. Worsley, J. P. Poline, C. D. Frith and R. S. J. Frackowiak, Statistical parametric maps in functional imaging: A general linear approach, Hum. Brain Mapp. 2(4) (1994) 189–210. 32. P. Piaggi, D. Menicucci, C. Gentili, G. Handjaras, A. Gemignani and A. Landi, Adaptive filtering and random variables coefficient for analyzing functional magnetic resonance imaging data, Int. J. Neural Syst. 23(3) (2013) 1350011. 33. D. Zhang and M. E. Raichle, Disease and the brain’s dark energy, Nat. Rev. Neurol. 6(1) (2010) 15–28.

1440008-15

A predictive modeling approach to analyze data in EEG-fMRI experiments.

In this paper, a novel technique based on blind source extraction (BSE) using linear prediction is proposed to extract rolandic beta rhythm from elect...
4MB Sizes 2 Downloads 5 Views