How to validate similarity in linear transform models of event-related potentials between experimental conditions?

Journal of Neuroscience Methods 236 (2014) 76–85

Contents lists available at ScienceDirect

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

Computational Neuroscience

How to validate similarity in linear transform models of event-related potentials between experimental conditions? Fengyu Cong a,b,∗ , Qiu-Hua Lin c , Piia Astikainen d , Tapani Ristaniemi b a

Department of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, China Department of Mathematical Information Technology, University of Jyväskylä, Finland School of Information and Communication Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, China d Department of Psychology, University of Jyväskylä, Finland b c

h i g h l i g h t s • Linear transform models (LTMs) of ERP data are investigated. • Using ICA, relative mapping coefficients (RMC) in LTMs are defined. • Using RMCs of an ERP, similarity in LTMs between conditions is examined.

a r t i c l e

i n f o

Article history: Received 14 May 2014 Received in revised form 18 August 2014 Accepted 18 August 2014 Available online 23 August 2014 Keywords: Event-related potentials Independent component analysis Linear transform model Mapping coefficient

a b s t r a c t Background: It is well-known that data of event-related potentials (ERPs) conform to the linear transform model (LTM). For group-level ERP data processing using principal/independent component analysis (PCA/ICA), ERP data of different experimental conditions and different participants are often concatenated. It is theoretically assumed that different experimental conditions and different participants possess the same LTM. However, how to validate the assumption has been seldom reported in terms of signal processing methods. New method: When ICA decomposition is globally optimized for ERP data of one stimulus, we gain the ratio between two coefficients mapping a source in brain to two points along the scalp. Based on such a ratio, we defined a relative mapping coefficient (RMC). If RMCs between two conditions for an ERP are not significantly different in practice, mapping coefficients of this ERP between the two conditions are statistically identical. Results: We examined whether the same LTM of ERP data could be applied for two different stimulus types of fearful and happy facial expressions. They were used in an ignore oddball paradigm in adult human participants. We found no significant difference in LTMs (based on ICASSO) of N170 responses to the fearful and the happy faces in terms of RMCs of N170. Comparison with existing method(s): We found no methods for straightforward comparison. Conclusions: The proposed RMC in light of ICA decomposition is an effective approach for validating the similarity of LTMs of ERPs between experimental conditions. This is very fundamental to apply grouplevel PCA/ICA to process ERP data. © 2014 Elsevier B.V. All rights reserved.

1. Introduction

∗ Corresponding author at: Department of Biomedical Engineering, Faculty of Electronic Information and Electrical Engineering, Dalian University of Technology, ChuangXinYuan Dasha, Room A0934, LingGong Road #2, Dalian 116024, China. Tel.: +86 158 411 92277. E-mail addresses: [email protected], [email protected] (F. Cong). http://dx.doi.org/10.1016/j.jneumeth.2014.08.018 0165-0270/© 2014 Elsevier B.V. All rights reserved.

EEG data are modeled as the linear transform model (LTM) under the EEG frequency range (Makeig et al., 1996, 1997, 1999). In this model, the data collected along the scalp are the mixtures of sources which are of electrical brain activity, and a mixing/mapping matrix connects the sources and the mixtures. The coefficients of the matrix contain mapping coefficients between sources in brain and points along the scalp. In an experiment to elicit event-related

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

potentials (ERPs), EEG data processing methods are often applied to estimate ERPs’ components in light of the LTM model (Luck, 2005; Sanei and Chambers, 2007). As an advanced signal processing method, independent component analysis (ICA) conforming to a LTM has been successfully applied to separate scalp EEG data into sources of ERPs (Makeig et al., 1996; Vigario and Oja, 2008). It has been performed on the EEG data of a single trial (Cong et al., 2010; Iyer and Zouridakis, 2007), the concatenated EEG data of a number of single trials (Delorme and Makeig, 2004; Eichele et al., 2011), and the averaged EEG data over many single trials (Cong et al., 2011b; Makeig et al., 1997). In this study, we focus on the data processing for the averaged EEG data. For example, in a passive oddball paradigm to elicit ERPs using pictures of human faces as stimuli, the happy and the fearful expressions have been applied as infrequently presented ‘deviant’ stimuli among natural ‘standard’ faces (Astikainen et al., 2013; Astikainen and Hietanen, 2009). Then, for a participant’s averaged EEG data, the LTM regarding an electrode site for two deviants reads xm,h,l (t) = am,1,h,l S1,h,l (t) + · · · + am,n,h,l Sn,h,l (t) + · · · + am,N(l),h,l SN(l),h,l (t) + vm,h,l (t)

xm,f,l (t) = am,1,f,l S1,f,l (t) + · · · + am,n,f,l Sn,f,l (t) + · · · + am,N(l),f,l SN(l),f,l (t) + vm,h,l (t) where ‘h’ denotes the deviant of the happy expression, and ‘f’ represents the deviant of the fearful expression, ‘l’ is the index for a participant, ‘n’ symbolizes the number of the source, N(l) is the number of all sources for the participant #l, and ‘m’ is the number of the electrode. For the lth participant under the deviant of the happy expression, xm,h,l (t) denotes the averaged EEG data, Sn,h,l (t) represents the nth source of electrical brain activity, and am,n,h,l is the mapping coefficient for Sn,h,l (t) to the point where the electrode #m is placed, and m ∈ [1, M], n ∈ [1, N(l)], n ∈ [1, N(l)], and l ∈ [1, L], M is the number of all electrodes, and L is the number of all participants. It should be noted that am,n,h,l depends on the properties of the volume conductor in brain of the participant #l, the location of the nth source in brain, and the measurement point along the scalp (Makeig et al., 1999). For the averaged EEG data, group ICA is often applied with the concatenation of ERP data of different experimental conditions (Kalyakin et al., 2009, 2008) or/and different participants (Kovacevic and McIntosh, 2007; Vakorin et al., 2010). Using group ICA, only one set of independent components and an unmixing matrix are estimated (Eichele et al., 2011). Such an approach does inherently assume that the LTMs of different experimental conditions or/and different participants are theoretical identical (Cong et al., 2013b). This means the mapping coefficients of a source keep identical across different experimental conditions or/and participants, i.e., am,n,h,l = am,n,f,l , or/and am,n,h,l1 = am,n,h,l2 , ∀li ∈ [1, L], i = 1, 2, l1 = / l2 , and the orders of sources remain the same along different experimental conditions or/and participants (Cong et al., 2013b). They are indeed the general hypotheses in the electrical fields of the brain (Nunez and Srinivasan, 2005). However, within the knowledge of the present authors, we have not found any previous study to straightforwardly validate these hypotheses in terms of signal processing methods despite that ICA (Eichele et al., 2011), principal component analysis (PCA) (Dien, 2012; Lohvansuu et al., 2013) and tensor decomposition (Cong et al., 2013c) have utilized such assumption. This motivates us to formulate a signal processing approach for validating the hypotheses. In this study, we develop an ICA-based procedure to examine the similarity in LTMs of an ERP, N170, between two stimulus types (fearful faces vs happy faces) in the example mentioned above.

77

2. Method 2.1. Data description Twenty two healthy adult volunteers were recruited by a newspaper advertisement (18 females, age range 30–58 years, mean 46.1 years) to take part in the experiment for data collection. All participants were right-handed and reported to have normal or corrected-to-normal vision. An informed written consent was obtained from each participant. The experiment was undertaken in accordance with the Declaration of Helsinki, and the ethical committee of the University of Jyväskylä approved the research protocol. In order to elicit ERPs, pictures of emotionally expressive faces were presented to the participants. The stimuli were neutral, fearful and happy faces of four different actors from the series ‘Pictures of Facial Affect’ (Ekman and Friesen, 1976). The stimuli were presented in a passive oddball condition: the pictures of neutral facial expressions served as a repeated standard stimulus (probability = 0.8) and the pictures of the happy and fearful expressions (probability = 0.1 for each) as rarely presented deviant stimuli. During the recordings, the participants sat in a chair, and were instructed to pay no attention to the visual stimuli but instead focused on a radio play presented via loud speakers. At least two standards were presented between randomly presented consecutive deviants. The stimulus duration was 200 ms, and the stimulus onset asynchrony was 700 ms. Such a paradigm elicited a so called N170 response which is sensitive for faces and enhanced for emotional faces compared to neutral faces (Astikainen et al., 2013; Astikainen and Hietanen, 2009). EEG data were collected with 14 electrodes including Fz, F3, F4, Cz, C3, C4, Pz, P3, P4, P7, P8, Oz, O1 and O2 according to the international 10–20 system through Brain Vision Recorder software (Brain Products GmbH, Munich, Germany). An average reference was used. Bipolar electrodes were placed above and below the left eye and lateral to the left and right orbit to measure the eye movements and blinks. Data were digitally on-line filtered from 0.1 to 100 Hz. The sampling frequency was 1000 Hz. Data were first offline preprocessed with Brain Vision Analyzer software (Brain Products GmbH, Munich, Germany). Continuous data were segmented (from 200 ms pre-stimulus period to 500 ms after the stimulus onset) and the baseline was corrected against 200 ms pre-stimulus interval. Segments with amplitude values beyond the range between −100 and 100 ␮V in any recording channel, including the electrooculogram channel, were rejected. The number of kept trials for the averaging was about 100 per deviant type. After the data preprocessing, the remaining trials were averaged to produce the ERP waveforms, i.e., the averaged EEG data. 2.2. Procedure of ICA to estimate peak amplitude of an ERP After the conventional ERP data processing, the peak amplitude of an ERP is often measured from the averaged EEG for the further statistical analysis. As shown in the Introduction, the averaged EEG data are still mixtures of many brain sources. Using ICA, the averaged EEG data can be spatially filtered to theoretically produce the sole waveform of one brain source in the electrode field (Cong et al., 2013a, 2011c, 2011d). 2.2.1. Complete ICA procedure as a spatial filter When sensor noise is omitted or included in the LTM (Cong et al., 2014), the LTM associating the EEG data (x) along the scalp and the electrical sources (s) in brain can be expressed as x = As

(1)

78

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

where x = [x1 , x2 , · · · , xM ]T , s = [s1 , s2 , · · · , sN ]T , and A ∈  M×N with the full column rank is usually called as the mixing matrix regarding ICA. In this study we name it as the mapping matrix. Each column of the matrix contains coefficients for mapping a source in brain to points along the scalp and the coefficients formulate the scalp topography of the source. For any source, its mapping can be illustrated as xr = ar · sr

(2) ]T ,

where xr = [x1r , x2r , · · · , xM,r ar is one column of A, and r ∈ [1,N]. In this case, xr is not the mixture like x in (1) any more, but only contain sole information of one brain source. Hence, one goal to apply ICA is to achieve the mapping of one source in (2) from the mixture in (1) (Cong et al., 2011c,d). For simplicity without losing generality, we assume M = N here. In order to obtain (2), an unmixing matrix is first learned by ICA (Comon, 1994; Hyvarinen et al., 2001), and then it transforms the mixture in (1) into independent components as y = Wx

(3) ]T

where y = [y1 , y2 , · · · , yN denotes all the extracted components by ICA. It should be noted that the variance and polarity of any component, yn , is not determined in (3). However, the peak amplitude of an ERP is very important parameter (Luck, 2005). Therefore, it is necessary to correct the indeterminacy when ICA is applied to study ERPs. After ICA decomposition, one component of interest is often selected according to prior knowledge of the desired brain activity. It is then projected back to the electrode field to correct the variance indeterminacy of the extracted component by ICA (Makeig et al., 1997, 1999; Onton et al., 2006) through (4)

B = W−1

(5) ]T ,

where ek = [e1,k , e2,k , · · · , eM,k bk is one column of B, yk is one element of y and k ∈ [1,N]. Consequently, we obtain the desired electrical brain activity’s determined magnitude with the unit of the microvolt in the context of EEG data (Makeig et al., 1997, 1999; Onton et al., 2006). Furthermore, with the global matrix of ICA (Cichocki and Amari, 2003), a component can be interpreted as N 

ckn sn

(6)

n=1

C = WA

(7)

where ckn is the (k,n) element of the global matrix C. Then, (4) turns to ek = bk ·

2.2.2.2. Over-determined model. From Eq. (3) to Eq. (10), the assumption is that the number of sensors is equal to the number of sources. When the high-dense array is used to collect EEG data, it tends to assume that the number of sensors is bigger than the number of sources. In this case, the model is over-determined, i.e., M > N. The dimension reduction is usually applied before ICA decomposition. This process reads y = WVx = WVAs

(11)

where V ∈  N×M

is learned by PCA and model order selection (Cong et al., 2011d, 2014). Then, the projection of one ICA component to the electrode field turns to

ek = tk · yk

(12)

T

tk = V · bk

(13)

where tk denotes the topography of the component #k. Therefore, under the theoretical global optimization, the projection is BC=A

ek = bk · yk

yk =

where ck,j is the nonzero element in the kth row and jth column of the global matrix C, aj is one column of A, xj = [x1,j , x2,j , · · · , xN,j ]T , and j ∈ [1,N]. Consequently, when ICA is globally optimized to separate EEG data, we still cannot obtain the true source or the real mixing matrix, but can theoretically gain the individual sources of brain activities with unknown scales as shown by (9) and the determined mapping of individual sources along the scalp as illustrated by (10). For the same source, (2) and (10) are identical. Eq. (9) has been generally well-known in the field of ICA and (10) has been seldom noticed. From (1) to (10), the EEG data are spatially filtered by the linear transform.



N n=1



ckn sn

(9)

Subsequently, (8) becomes (Cong et al., 2013a, 2011c,d) BC=A

ek = bk · (ckj sj ) = aj · sj = xj

(14)

2.2.3. Back-projection of one ICA component: practical expectation In contrast to (8) and (10), what we can obtain in practice just approximate the counterparts in theory since ICA is often locally optimized in fact (Himberg et al., 2004). This means that there are probably at least two or more nonzero elements in some columns or rows of global matrix C of ICA. Moreover, we cannot know the error between what we practically obtain and what we theoretically gain. Two principles should be acknowledged: (1) if multiple runs of ICA decomposition on EEG data are not stable, results are not acceptable (Himberg et al., 2004); (2) although we cannot know the true source of electrical brain activity, we may know the function of the brain activity; therefore, the extracted source by ICA from EEG data should be evaluated by the function of the corresponding brain activity to validate whether what has been extracted are reliable or not (Vigario and Oja, 2008).

(8)

2.2.2. Back-projection of one ICA component: theoretical expectation 2.2.2.1. Determined model. Here, we define that under the global optimization of ICA there is only one non-zero element in each row and each column in C and each component is just the scaled version of a source with the unknown scale (Cichocki and Amari, 2003; Hyvarinen et al., 2001). With such a condition taken into account, (6) reads yk = ckj sj

ek = VT · bk · (ckj sj ) = aj · sj = xj

(10)

2.2.4. Relative mapping coefficient (RMC) Now, we suppose that the ICA decomposition is globally optimized for averaged EEG data of each participant under each experimental condition to elicit N170 in this study. In terms of the projection of the extracted N170 component for the fear and the happy deviant for one participant, the projection can be expressed as em,f = am,f · sf

(15)

em,h = am,h · sh

(16)

where m = 1,. . .,M, sf , and sh are sources of N170 under the fear and the happy deviant, am,f and am,h are their mapping coefficients to the electrode #m. As mentioned earlier, it is generally assumed that am,f and am,h are identical in the electrical fields of the brain in

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

an ERP experiment (Nunez and Srinivasan, 2005). We attempt to validate such an assumption through a signal processing approach here. As interpreted above, we cannot obtain the mapping coefficients am,f and am,h using ICA even though ICA decomposition is globally optimized. However, we have found that the ratio of two mapping coefficients of an ERP between any two electrode sites (denoted by m1 and m2 , m1 ∈ [1, M], m2 ∈ [1, M]) under one experiment condition is available. Such a ratio can be then defined as am1 ,f am1 ,f · sf em1 ,f m1 m = = = (17) 2 ,f am2 ,f am2 ,f · sf em2 ,f m1 m = ,h 2

am1 ,h am2 ,h

am1 ,h · sh

=

am2 ,h · sh

=

em1 ,h

(18)

em2 ,h

m1 and  m1 must be idenIf ami ,f and ami ,h (i = 1, 2) are equal m ,f m ,h 2

2

m1 and  m1 tical too. Provided that m are the same, we cannot ,f m ,h 2

2

directly make the conclusion that ami ,f and ami ,f (i = 1, 2) are alike with the probability 100%. However, since the fear deviant and the happy deviant randomly appeared in the experiment, it is reasonable to expect that the two mapping coefficients of an ERP under two experiment conditions at the same electrode should be in the same quantitative level, i.e., the level of signal amplitude. Based on m1 and such hypothesis, we can come to the conclusion that if m ,f 2

m1 are identical, a m mi ,f and ami ,h (i = 1, 2) are almost equal. ,h 2

Hence, in order to reliably reveal such a ratio of two mapping coefficients of an ERP between any two electrode sites under one experiment condition, we define the relative mapping coefficient (RMC) as the following, m2 ,f =

m2 ,h =

1 M 1 M

M  

m1 m ,f



2

(19)

m1 =1 M  

m1 m ,h 2



(20)

m1 =1

where · is the absolute value of a variable, and m2 is the number of the electrode site, m2 ∈ [1, M]. Here, the RMC of an ERP under one condition at one electrode site here is the average over all the ratios between the mapping coefficient of this ERP at this electrode site and its mapping coefficients at all electrode sites. We calculate the absolute value of the ratio here because the polarities of an ERP in different regions of the scalp tend to be different. In this study, we compared RMC between two experimental conditions, i.e., the two deviant stimuli, given one electrode site. If the difference in RMCs of an ERP between two experimental conditions were not significant, it could be assumed that the mapping coefficients of this ERP under the two experimental conditions were very similar. 2.2.5. Six steps for data processing In this study, the data processing included six steps which are illustrated as the following. Step 1: Averaging. As mentioned earlier in the Method section, the paradigm in this study was the passively oddball paradigm. For such a paradigm, the signal to noise ratio is relatively low in contrast to the active oddball paradigm to elicit, for example, P300. In this case, it is better to perform ICA on the averaged trace, not the EEG data of concatenated single trials (Cong et al., 2011b). Hence, after the artifacts rejection, all kept trials were averaged to produce the averaged trace at each electrode site and each deviant type separately for each participant. Step 2: Single-channel filtering. In our previous study, we illustrated the benefits of an optimal single-channel filter, namely the optimal wavelet filter for the desired ERP (Astikainen et al., 2013;

79

Cong et al., 2012, 2011e). Due to only 14 sensors used in our study, the appropriately designed wavelet filter assisted to convert the absolutely underdetermined linear transform model of the ordinary averaged data to the quasi-determined one (Astikainen et al., 2013; Cong et al., 2011e). This is because the used wavelet filter was in light of the properties of the desired ERP and could remove many sources of no interest (Astikainen et al., 2013; Cong et al., 2011b). Hence, it was used here too. For the appropriately designed wavelet filter, the reversal biorthogonal wavelet with the order of 6.8 was used, the averaged EEG data were decomposed into ten levels, and the coefficients from the sixth and the ninth were used to reconstruct the desired ERP component (Astikainen et al., 2013; Cong et al., 2011e). Consequently, frequency responses of such a wavelet filter conform to spectral properties of the ERP elicited by the passively oddball paradigm (Astikainen et al., 2013), which is fundamental to wavelet based analysis for ERPs (Basar et al., 2001). Fig. 2 shows the frequency responses of the wavelet filter. Based on the left figure for magnitude responses, we can make three conclusions: (1) frequency contents from about 1.5 to 10 Hz in the original ERP waveform are almost entirely kept in the filtered ERP waveform, (2) frequency contents from about 1.5 to 0 Hz in the original ERP waveform are decreased from 0 dB to about 20 dB, (3) frequency contents from 10 to 22 Hz in the original ERP waveform are gradually decreased from 0 dB to about 20 dB. According to the right figure, the linear phase responses mean that the ERP waveform cannot be severely distorted by the filter. Hence, the wavelet filter used in our study can at least reduce some of the fluctuating brain activities (for example, alpha, beta, etc.) from the ERP waveform, as well as the sensor noise. By this way, some brain activities can be removed, and then, the number of brain activities can be reduced. This is very important and fundamental for the further ICA application (Cong et al., 2011e). Step 3: Multi-channel decomposition. ICASSO (Himberg et al., 2004) has been proven to be an effective approach for implementing ICA decomposition in extracting ERPs elicited by the passively oddball paradigm (Cong et al., 2011b; Kalyakin et al., 2009, 2008). It runs one ICA algorithm many times respectively with individually and randomly initialized unmixing matrices. Then, all the extracted components are clustered into the predefined number of clusters. Finally, each common component in each cluster represents one component extracted by ICASSO and the stability index denoted by Iq is calculated for such a component (Himberg et al., 2004). The Iq ranges from ‘0’ to ‘1’. When approaching to ‘1’, it indicates that the corresponding component is extracted out in almost every round of ICA decomposition. This implies the high stability of the ICA decomposition and satisfactory reliability of that component. Otherwise, it means the ICA decomposition is not stable and then the estimated component cannot be reliable. The magnitude of Iq has been regarded as the criterion to evaluate the ICA decomposition on brain signals (Cong et al., 2011b,e). For ICASSO, the FastICA algorithm (Hyvarinen, 1999) was run 100 times and 14 components were extracted at each run, and the default sets for other parameters and algorithms were used (Himberg et al., 2004). It should be noted that ICASSO was performed on the data under each deviant for each participant individually in this study. Due to 14 electrodes, we extracted 14 components under each run of ICA decomposition in ICASSO. Since the determined ICA algorithm was used we indeed assumed 14 sources in the filtered data. If the ICA decomposition was stable, such assumption should be reasonable, as well as acceptable in practice. Step 4: Desired ERP component selection. As illustrated by (4), one component of interest may be selected for the further processing. In this study, we selected to study N170 component. Hence, after ICA decomposition, we chose the component possessing the sole peak with the latency around 170 ms as the desired one.

80

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

Fig. 1. Frequency responses of the wavelet filter used in this study.

For each of all the 22 participants under each deviant, there always existed such a component among all the 14 components extracted by ICA. Step 5: Projection of the selected component. As explained by (4), the selected component is projected back to the electrode field. Step 6: Analyzing relative mapping coefficient: After the back projection of N170 component in the electrode field is obtained, the RMC can be calculated according to (19) and (20) under each deviant at each electrode for each participant. And then, through ANOVA (analysis of variance), we examined the difference of RMCs between two deviants. 3. Results Here, we first demonstrate the reliability of the ICA decomposition which is the basis for the further analysis, and then, we present the similarity in RMCs between two deviants to elicit N170 in this study. 3.1. Reliability of ICA decomposition Fig. 2 shows the stability (denoted by Iq) of ICA decomposition for both deviants of each participant. Indeed, each Iq in Fig. 1 is the averaged Iq over 14 components for each deviant and participant. What we are interested is just N170 here. However, the non-reliable estimation of other components also affects the projection of the desired component extracted by ICA in the electrode field (Cong et al., 2011c,d). Hence, we used the averaged Iq over 14 components to represent the stability of ICA decomposition in our study. Fig. 2 reveals that such an Iq is over 0.85 in most cases, and this means the ICA decomposition is stable. As mentioned earlier, except the performance of ICA decomposition, it is also necessary to analyze how plausible the estimated ERP by ICA is in light of the psychophysiological interpretation (Vigario and Oja, 2008). Topography of an ERP is one of the most important characteristics to ascertain this issue (Blankertz et al., 2011; Delorme and Makeig, 2004). Regarding N170 of healthy adults, it is observed maximally at lateral parietal sites. The N170

Table 1 Mean amplitudes of the back-projected ICA components and statistical tests comparing the difference of the peak amplitude of N170 between left and right hemisphere in the parieto-occipital area. Method|Mean amplitude (␮V)

Left

Right

F(1,21)

p

Ordinary average-fear deviant Ordinary average-happy deviant ICA-fear deviant ICA-happy deviant

−2.4850 −2.5416 −1.8128 −1.6462

−3.3046 −3.6040 −2.4138 −2.4164

11.3570 11.0318 12.1177 10.3665

0.0029 0.0032 0.0022 0.0041

amplitude is typically larger in the right than the left parietooccipital areas. Now, N170 amplitude estimated by ICA should possess such a property. Fig. 3 presents the ERPs achieved by the conventional grand averaging of the EEG data. Fig. 4 shows the grand averages of N170 estimated by ICA. For the statistical tests, the peak amplitudes of N170 at P7 and O1 in each participant were averaged to represent the amplitude of N170 at the left parieto-occipital area. Correspondingly, the N170 at P8 and O2 were averaged to represent the right parieto-occipital response. Responses to happy and fearful faces were analyzed separately. As shown in Table 1, the N170 magnitude was significantly larger in the right than left parieto-occipital channels both in the averages based on the raw data as well as in those estimated by ICA. 3.2. Analysis of peak amplitude of N170 The grand averaged peak amplitude of N170 is shown in Figs. 3 and 4. This is very conventional to show in an ERP study (Luck, 2005). In order to examine the difference in the peak amplitudes of N170 between happy and fearful faces, the amplitudes at P7, P8, O1 and O2 were statistically analyzed. These four sites are typical for N170 (Rossion and Jacques, 2008). Table 2 presents the statistical tests of those amplitudes. For the amplitudes measured from the conventionally preprocessed waveforms, the difference between two emotional faces was significant at O2 and almost at O1. For the amplitudes measured from the waveforms processed by the ICA approach, the difference disappeared. In contrast to O1 and O2, P7 and P8 are more typical

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

Fig. 2. Stability (denoted by the averaged Iq over all 14 extracted components) of ICA decomposition for each participant under each deviant.

Fig. 3. Grand ordinary average of preprocessed EEG recordings. The preprocessing is introduced in the third paragraph of Section 2.1.

81

82

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

Fig. 4. Grant average of N170 estimated by ICA in the electrode field.

Table 2 Statistical tests of difference in the peak amplitudes of N170 between two emotional faces at typical electrode sites. Parameter|electrode

P7

P8

O1

Statistical tests for amplitudes of N170 by conventional method 0.9 1.6 5.3 F(1,21) 0.352 0.224 0.052 p Statistical tests for amplitudes of N170 by the ICA approach F(1,21) 2.090 0.640 0.170 p 0.163 0.433 0.684

O2 11.3 0.003 0.534 0.473

sites for N170 (Rossion and Jacques, 2008). No matter why method was applied, the difference in the amplitudes of N170 was not significant between two emotional faces at P7 and P8. For an ERP, the effects between experimental conditions can be more reliably observed from the data at more typical sites (Rossion and Jacques, 2008). Therefore, we believe that in our study the two emotional faces elicited N170 without significant different amplitudes. The reason for the significant difference in the amplitudes of N170 at O1 and O2 in the conventionally processed waveforms is that some interference might exist at those two sites. After ICA was applied, the real N170 source is extracted and the difference at O1 and O2 therefore disappeared in our study.

3.3. Analysis of relative mapping coefficient Fig. 5 presents the RMC of N170 in each participant at each electrode site for both the happy and the fearful faces. For most participants at most electrode sites, the RMCs for the happy and the fearful faces were very similar. Statistically, Table 2 demonstrates that the difference of the RMCs between the two facial expressions was not significant at any electrode site. Such results well correspond to the analysis of the peak amplitudes of N170 estimated by the ICA approach (Table 3).

4. Discussion In this study, we proposed one ICA-based approach to validate the similarity in linear transform models of an ERP between two experiment conditions. With this approach, we found that the linear transform models of visual N170 of the healthy adults between the fear and the happy deviant in a passively oddball paradigm were similar. This does keep consistent to the general assumption in the electrical fields of the brain in an ERP experiment (Nunez and Srinivasan, 2005). ICA has become an important tool to analyze brain signals (Delorme and Makeig, 2004; Vigario and Oja, 2008). Regarding group ICA, data of responses from different experiment stimuli, for example, different deviants in an oddball paradigm, are often concatenated (Kalyakin et al., 2009, 2008). Such concatenation theoretically assumes that the electrical fields of the brain keep identical among different external stimuli in an experiment. However, before our ICA-based approach, we find no reports to straightforwardly validate this assumption in terms of signal processing. Briefly to say, the application of ICA to study ERPs includes two steps. The first is to decompose the ERP data, i.e., separate the mixtures, for obtaining the independent components. Secondly, projecting the extracted component of interest back into the electrode field recovers the true scale of EEG with the unit of microvolt for comparison (Makeig et al., 1997, 1999). We have found that theoretically what we can obtain is neither the true source, nor the mapping coefficient of the true source from the location in the brain to the point along the scalp, but the multiplication of them (Cong et al., 2011c,d). Consequently, we cannot directly compare the mapping coefficients representing the linear transform models of an ERP between two experiment conditions, for example, two deviants in one oddball paradigm. Furthermore, most of research using ICA to study an ERP did not validate whether the projection really met the theoretical expectation of this ERP in the electrode field or not (Delorme and Makeig, 2004; Kalyakin et al.,

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

83

Fig. 5. Relative mapping coefficient of N170 for each participant under each deviant.

2009, 2008; Marco-Pallares et al., 2005; Vigario and Oja, 2008). Our previous studies showed that the polarities of the projection of one extracted independent component in the electrode field may not be fully corrected due to the inevitable local optimization of ICA decomposition in practice (Cong et al., 2011b,c,d). This means when yk is closely associated with the desired brain activity that the source sj corresponds to, the sign of ei,k in (4) may be probably different with that of xi,j in (10) (i ∈ [1,N]). Therefore, this is one criterion to examine whether the ICA decomposition is satisfactory or not for an ERP component. If the polarity of a well-known ERP in the back-projection in the electrode field was reversed at a typical point along the scalp, the ICA decomposition for that ERP

component would not be accepted. We examined all the backprojection of N170 components of all participants in the electrode field at P7 and P8 (typical electrode locations for N170 (Rossion and Jacques, 2008), and found that no polarity reversal occurred in the back-projection. It should be noted that the RMC defined here can be used in the determined LTM and the over-determined model. In this study, the 14-senosor-14-source model was applied in this study since only 14 sensors were used for data collection. Such a number is comparable to some previous studies. For example, we found that 11 sources were estimated by the BIC method in an ERP experiment when group ICA was applied (Kovacevic and McIntosh, 2007). In our

Table 3 Mean values of relative mapping coefficient of N170 and statistical tests of difference in relative mapping coefficient of N170 between two deviants at each electrode. Parameter|electrode

Fz

Oz

F3

F4

C3

C4

P3

RMC-happy RMC-fear F(1,21) p

0.5146 0.2329 1.8080 0.1931

0.3039 0.2503 0.2446 0.6260

0.2083 0.4194 1.7432 0.2009

0.1838 0.3933 1.4012 0.2497

0.3981 0.1504 2.0309 0.1688

0.3342 0.1201 1.4207 0.2466

0.1622 0.0848 0.9547 0.3396

Parameter|electrode

P4

Cz

Pz

P7

P8

O1

O2

RMC-happy RMC-fear F(1,21) p

0.0962 0.0728 2.3383 0.1412

0.5907 0.2543 1.7709 0.1975

0.3121 0.1602 1.6285 0.2158

0.5783 0.3468 1.6375 0.2146

0.9574 0.3947 2.0878 0.1632

0.3124 0.1362 1.8205 0.1916

0.4119 0.2230 1.5673 0.2244

84

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85

previous studies, about 14 sources were estimated by the model order selection method called GAP in two ERP experiments after the ERP data were filtered by the wavelet filters (Cong et al., 2011a, 2013b). Regarding the electrode fields of the brain, we assume that the mapping coefficients of an ERP between two deviants in one oddball paradigm are in the same quantitative level at one electrode site. And then, we defined the relative mapping coefficient of an ERP under one experiment condition, i.e., one deviant in our study, as shown by (17)–(20). If the difference of RMCs of an ERP between two deviants is not significant in practice, it is reasonable to predict the mapping coefficients of this ERP between two deviants are very probably identical in theory. The estimation of RMC in this study was reliable since we validated that the estimation of N170 was reliable both in terms of ICA decomposition and in light of psychophysiology properties of N170 as stated in the Results. Indeed, RMC can be very susceptive to un-reliable ICA decomposition since it is the ratio between the projections of one component at two electrode sites. Any component which is not reliably extracted by ICA may be a big risk of affecting the projection of any other components in the electrode field. The negative effect of such situation is that the magnitude of the projection may become exceptionally large or small at some electrode sites (Cong et al., 2011c,d), resulting in unusually large or small RMC consequently. Hence, before measuring RMC, it is necessary to validate the effectiveness of the estimation by ICA. Otherwise, RMC is not acceptable. Furthermore, what we can validate between two experiment conditions in this study is just one column of the mixing matrix in (1) because we only analyzed one ERP, N170 in this study. Indeed, in one ERP experiment it can be very challenging to analyze all ERP components extracted by ICA due to the following two reasons. Firstly, an ERP experiment is usually designed to elicit one particular ERP component and other ERPs inevitably elicited by the experiment are not as stable as the prime one. This does bring difficulty for ICA to as well extract those ERP components as the prime one from every participant individually. Secondly, it is impossible that all the ERPs elicited in one ERP experiment can be well illustrated. In order to validate the effectiveness of the estimation of an ERP through ICA, it is necessary to interpret the estimation based on the psychophysiology knowledge of this ERP (Vigario and Oja, 2008). We cannot guarantee that all the ERP components extracted by ICA are reliable from this point of view despite that the estimation of those components is stable in terms of ICA decomposition. An interesting early study developed an effective method to determine the similarity of the underlying linear mixing processes of multiple participants/single trials (Hyvarinen, 2011). For example, in this study, we extracted 14 temporal components from the 14-sensor EEG data for each participant for each deviant. Each component is associated with one scalp topography which is denoted by one column of the mixing matrix in Eq. (1) in the manuscript. Since 22 adults participated in the experiment, for two deviants we could altogether obtain 616 (14 by 2 by 22) topographies representing by a matrix with the sizes of 616 by 14. After the matrix is normalized columnwisely, we could perform a clustering method on the matrix to cluster the 616 topographies into the pre-defined number of clusters. Ideally, in each cluster, the topographies belonging to one cluster can be very similar with each other. Because the label of each scalp topography is known, we are able to obtain how many participants share similar topographies for both deviants. However, since each scalp topography is associated with one temporal component in ICA, it is hard to guarantee that the temporal components associated with the topographies within one cluster all correspond to the same ERP component, say N170 in this study. For RMC, the temporal components should belong to the same ERP component. This is the key difference between the proposed RMC and the approach

developed by Hyvarinen (2011). Hence, we do not apply that approach in our study. Moreover, there are three open questions which are worth further investigation in the future. The first is that we cannot determine which column of the mixing matrix in (1) can be examined by the proposed RMC due to the permutation indeterminacy of ICA. The order of the extracted components by ICA is not determinate inherently (Hyvarinen et al., 2001). Secondly, we cannot determine whether the position of the column corresponding to the ERP source in the mixing matrix under one experiment condition is the same to that under another condition. At last, in this study the subject was the participant in the experiment. If we estimated N170 from each single trial, we could even validate the similarity of linear transform models between two experiment conditions of an ERP in the level of single trials with the proposed RMC. Indeed, it will be very challenging to estimate reliable ERPs by ICA from each single-trial’s EEG data elicited by the passively oddball paradigm. We will try to answer those questions in the future publication. Acknowledgements Cong thanks the financial support from TEKES in Finland (MUSCLES projects, 40334/10), the Fundamental Research Funds for the Central Universities [DUT14RC(3)037] in Dalian University of Technology in China, and National Natural Science Foundation of China (Grant No. 81471742). All authors appreciate the anonymous reviewers for the invaluable comments to improve the study. References Astikainen P, Cong F, Ristaniemi T, Hietanen JK. Event-related potentials to unattended changes in facial expressions: detection of regularity violations or encoding of emotions? Front Hum Neurosci 2013;7:557, http://dx.doi.org/10.3389/fnhum.2013.00557; 10.3389/fnhum.2013.00557. Astikainen P, Hietanen JK. Event-related potentials to task-irrelevant changes in facial expressions. Behav Brain Funct 2009;5:30, http://dx.doi.org/10.1186/1744-9081-5 30. Basar E, Schurmann M, Demiralp T, Basar-Eroglu C, Ademoglu A. Event-related oscillations are ‘real brain responses’ – wavelet analysis and new strategies. Int J Psychophysiol: Off J Int Organizat Psychophysiol 2001;39(2–3):91–127. Blankertz B, Lemm S, Treder M, Haufe S, Muller KR. Single-trial analysis and classification of ERP components – a tutorial. Neuroimage 2011;56(2):814–25, http://dx.doi.org/10.1016/j.neuroimage.2010.06.048. Cichocki A, Amari S. Adaptive blind signal and image processing: learning algorithms and applications (Vol. revised). Chichester: John Wile & Sons Inc; 2003. Comon P. Independent component analysis, a new concept? Signal Process 1994;36(3):287–314. Cong F, Alluri V, Nandi AK, Toiviainen P, Fa R, Abu-Jamous B, Ristaniemi T. Linking brain responses to naturalistic music through analysis of ongoing EEG and stimulus features. IEEE Trans Multimed 2013a;15(5):1060–9. Cong F, He Z, Hamalainen J, Cichocki A, Ristaniemi T. Determining the number of sources in high-density EEG recordings of event-related potentials by model order selection. In: Proceedings on IEEE Workshop on Machine Learning for Signal Processing (MLSP) 2011; 2011a. p. 1–6. Cong F, He Z, Hamalainen J, Leppanen PHT, Lyytinen H, Cichocki A, Ristaniemi T. Validating rationale of group-level component analysis based on estimating number of sources in EEG through model order selection. J Neurosci Methods 2013b;212(1):165–72. Cong F, Huang Y, Kalyakin I, Li H, Huttunen-Scott T, Lyytinen H, Ristaniemi T. Frequency response based wavelet decomposition to extract children’s mismatch negativity elicited by uninterrupted sound. J Med Biol Eng 2012;32(3):205–14. Cong F, Kalyakin I, Ahuttunen-Scott T, Li H, Lyytinen H, Ristaniemi T. Single-trial based independent component analysis on mismatch negativity in children. Int J Neural Syst 2010;20(4):279–92. Cong F, Kalyakin I, Li H, Huttunen-Scott T, Huang YX, Lyytinen H, Ristaniemi T. Answering six questions in extracting childrenÂD ÂTM s mismatch negativity through combining wavelet decomposition and independent component analysis. Cognitive Neurodyn 2011b;5(4):343–59. Cong F, Kalyakin I, Ristaniemi T. Can back-projection fully resolve polarity indeterminacy of ICA in study of ERP? Biomed Signal Process Control 2011c;6(4):422–6. Cong F, Kalyakin I, Zheng C, Ristaniemi T. Analysis on subtracting projection of extracted independent components from EEG recordings. Biomed Technik/Biomed Eng 2011d;56(4):223–34. Cong F, Leppanen PH, Astikainen P, Hamalainen J, Hietanen JK, Ristaniemi T. Dimension reduction: additional benefit of an optimal filter for independent

F. Cong et al. / Journal of Neuroscience Methods 236 (2014) 76–85 component analysis to extract event-related potentials. J Neurosci Methods 2011e;201(1):269–80, http://dx.doi.org/10.1016/j.jneumeth.2011.07.015. Cong F, Phan AH, Astikainen P, Zhao Q, Wu Q, Hietanen JK, Cichocki A. Multi-domain feature extraction for small event-related potentials through nonnegative multi-way array decomposition from low dense array EEG. Int J Neural Syst 2013c;23:1–18, http://dx.doi.org/10.1142/S0129065713500068 (2(1350006)). Cong F, Puolivali T, Alluri V, Sipola T, Burunat I, Toiviainen P, Ristaniemi T. Key issues in decomposing fMRI during naturalistic and continuous music experience with independent component analysis. J Neurosci Methods 2014;223:74–84, http://dx.doi.org/10.1016/j.jneumeth.2013.11.025. Delorme A, Makeig S. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J Neurosci Methods 2004;134(1):9–21, http://dx.doi.org/10.1016/j.jneumeth.2003.10. 009. Dien J. Applying principal components analysis to event-related potentials: a tutorial. Dev Neuropsychol 2012;37(6):497–517. Eichele T, Rachakonda S, Brakedal B, Eikeland R, Calhoun VD. group EEGIFT: independent component analysis for eventEEG data. Comput Intell Neurosci 2011;2011:129365, related http://dx.doi.org/10.1155/2011/129365; 10.1155/2011/129365. Ekman P, Friesen WV. Pictures of facial affect. Palo Alto, CA: Consulting Psychologists Press; 1976. Himberg J, Hyvarinen A, Esposito F. Validating the independent components of neuroimaging time series via clustering and visualization. Neuroimage 2004;22(3):1214–22, http://dx.doi.org/10.1016/j.neuroimage.2004.03.027. Hyvarinen A. Fast and robust fixed-point algorithms for independent component analysis. IEEE Trans Neural Netw/Publ IEEE Neural Netw Council 1999;10(3):626–34, http://dx.doi.org/10.1109/72.761722. Hyvarinen A. Testing the ICA mixing matrix based on inter-subject inter-session consistency. Neuroimage 2011;58(1):122–36, or http://dx.doi.org/10.1016/j.neuroimage.2011.05.086. Hyvarinen A, Karhunen J, Oja E. Independent component analysis. New York: John Wile & Sons Inc; 2001. Iyer D, Zouridakis G. Single-trial evoked potential estimation: comparison between independent component analysis and wavelet denoising. Clin Neurophysiol: Off J Int Fed Clin Neurophysiol 2007;118(3):495–504, http://dx.doi.org/10.1016/j.clinph.2006.10.024. Kalyakin I, Gonzalez M, Ivannikov I, Lyytinen H. Extraction of the mismatch negativity elicited by sound duration decrements: a comparison of three procedures. Data Knowledge Eng 2009;68(12):1411–26.

85

Kalyakin I, Gonzalez N, Karkkainen T, Lyytinen H. Independent component analysis on the mismatch negativity in an uninterrupted sound paradigm. J Neurosci Methods 2008;174(2):301–12, http://dx.doi.org/10.1016/j.jneumeth.2008.07. 012. Kovacevic N, McIntosh AR. Groupwise independent component decomposition of EEG data and partial least square analysis. Neuroimage 2007;35(3):1103–12, http://dx.doi.org/10.1016/j.neuroimage.2007.01.016. Lohvansuu K, Hamalainen JA, Tanskanen A, Bartling J, Bruder J, Honbolygo F, Leppanen PHT. Separating mismatch negativity (MMN) response from auditory obligatory brain responses in school-aged children. Psychophysiology 2013. Luck SJ. An introduction to the event-related potential technique. The MIT Press; 2005. Makeig S, Bell AJ, Jung TP, Sejnowski TJ. Independent component analysis of electroencephalographic data. Adv Neural Inform Process Syst 1996;8:145–51. Makeig S, Jung TP, Bell AJ, Ghahremani D, Sejnowski TJ. Blind separation of auditory event-related brain responses into independent components. Proc Natl Acad Sci USA 1997;94(20):10979–84. Makeig S, Westerfield M, Jung TP, Covington J, Townsend J, Sejnowski TJ, Courchesne E. Functionally independent components of the late positive event-related potential during visual spatial attention. J Neurosci: Off J Soc Neurosci 1999;19(7):2665–80. Marco-Pallares J, Grau C, Ruffini G. Combined ICA-LORETA analysis of mismatch negativity. Neuroimage 2005;25(2):471–7, http://dx.doi.org/10.1016/j.neuroimage.2004.11.028. Nunez P, Srinivasan R. Electric fields of the brain: the neurophysics of EEG, 2. New York: Oxford University Press; 2005. Onton J, Westerfield M, Townsend J, Makeig S. Imaging human EEG dynamics using independent component analysis. Neurosci Biobehav Rev 2006;30(6):808–22, http://dx.doi.org/10.1016/j.neubiorev.2006.06.007. Rossion B, Jacques C. Does physical interstimulus variance account for early electrophysiological face sensitive responses in the human brain? Ten lessons on the N170. NeuroImage 2008;39(4):1959–79, http://dx.doi.org/10.1016/j.neuroimage.2007.10.011. Sanei S, Chambers JA. EEG signal processing. Wiley; 2007. Vakorin VA, Kovacevic N, McIntosh AR. Exploring transient transfer entropy based on a group-wise ICA decomposition of EEG data. Neuroimage 2010;49(2):1593–600, http://dx.doi.org/10.1016/j.neuroimage.2009.08.027. Vigario R, Oja E. BSS and ICA in neuroinformatics: from current practices to open challenges. IEEE Rev Biomed Eng 2008;1:50–61.