Cooperative Particle Filtering for Tracking ERP Subcomponents from Multichannel EEG

Abstract

In this study, we propose a novel method to investigate P300 variability over different trials. The method incorporates spatial correlation between EEG channels to form a cooperative coupled particle filtering method that tracks the P300 subcomponents, P3a and P3b, over trials. Using state space systems, the amplitude, latency, and width of each subcomponent are modeled as the main underlying parameters. With four electrodes, two coupled Rao-Blackwellised particle filter pairs are used to recursively estimate the system state over trials. A number of physiological constraints are also imposed to avoid generating invalid particles in the estimation process. Motivated by the bilateral symmetry of ERPs over the brain, the channels further share their estimates with their neighbors and combine the received information to obtain a more accurate and robust solution. The proposed algorithm is capable of estimating the P300 subcomponents in single trials and outperforms its non-cooperative counterpart.

1. Introduction

Event-related potentials (ERPs) are the electrical responses of the brain that directly measure the brain response to specific sensory, cognitive, or motor events. ERPs are typically generated as a result of a peripheral or external stimulation and are time locked to the stimulus. They are small voltage fluctuations (1–30 $\mathsf{\mu }$V) in the background EEG activity. ERPs can be used to assess several neurological disorders or cognitive processes. As the ERP components are quite small compared to the normal rhythm of the EEG signal, one common approach to extract them is by averaging time-locked single trials. In such an approach, it is assumed that the ERP components are constant over different trials and therefore, averaging over single trials attenuates the background EEG signal which behaves as random process here. However, this averaging scheme can result in loss of information related to trial-to-trial variability. This could be crucially important in several real scenarios where the ERP waveform changes over different trials due to habituation, fatigue, or the level of attention. Moreover, several brain abnormalities such as schizophrenia and depression can be identified by ERP variation over different trials [1,2]. Hence, it is important to develop an effective single trial estimation of the ERPs in order to consider the inter-trial variability of these waveforms. Moreover, accurate detection and estimation of single trial ERPs in time domain can have a significant impact on certain applications using entropy based analysis of ERPs which have not been evaluated extensively [3]. The use of multivariate multiscale entropy has great potential for dynamic analysis of ERPs in certain clinical applications such as in the discrimination of early and late stages of Alzheimer’s disease [4,5,6]. In addition, accurate detection of ERPs can lead to better separation of background EEG, which through entropy analysis, may reveal important information for applications such as dementia and schizophrenia.

Several methods have been proposed in the literature to analyze and estimate ERP components over single trials including Wiener [7], maximum a posteriori (MAP) [8], Kalman filtering [9,10], principal component analysis (PCA) [11], and independent component analysis (ICA) [12]. However, these methods might fail in the cases where the ERPs vary from trial to trial or the ERP signal to noise ratio (SNR) is low. Therefore, it is of interest to develop an effective and robust ERP analysis based on single trial estimation where the inter-trial variability is also taken into consideration.

One of the main components of ERP is the P300 wave which is a positive ERP component and occurs with a latency of about 300 ms after the stimulus. This component has been found to be associated with cognitive information processing involved in contextual updating, orientation of attention, and response resolution [13]. It has been observed that increasing the cognitive difficulty of the task results in a decrease of P300 amplitude and an increase in its latency [14,15,16]. The P300 wave consists of two main overlapping subcomponents with temporal correlations, known as P3a and P3b [16]. P3a has a frontocentral distribution since it is mainly generated from the prefrontal, frontal, and anterior temporal brain regions [17]. The P3a subcomponent is an unintentional response reflecting an automatic orientation of attention to a salient stimulus which is independent of the task in hand [17,18]. On the other hand, P3b is mainly biased toward the centroparietal region since it is generated principally by the posterior temporal, parietal, and posterior cingulate cortex regions [19]. Compared to P3b, P3a has a shorter latency and a more rapid habitation [17,19].

P300 subcomponents are often studied for the investigation of various psychiatric and neurological conditions, especially for the cases where the averaged P300 does not necessarily correspond to revelation of a particular condition such as mental fatigue [20,21]. Therefore, it is necessary to develop reliable and robust methods that not only use P300 amplitude and latency but also exploit the feature parameters related to P3a and P3b subcomponents. The method can then be reliably used for clinical diagnosis, assessment of mental activities, physical and mental fatigue, and in human-computer interface applications. In previous studies, blind source separation (BSS) and PCA have been applied to decompose the P300 waveform into its subcomponents [13,21,22,23,24]. Some other studies consider the averaged ERP and apply PCA to it, which is suitable for stationary data [25,26,27]. However, these methods might perform poorly or only estimate one of the components in the cases where there is a high temporal correlation between the subcomponents, correlated noise, or low ERP SNR.

In this study, we propose a robust method to separate and track the P300 subcomponents in single trial recordings. We benefit from the well-known symmetric distribution of the P300 waveform over the left and right hemispheres of the brain [28,29,30,31], and introduce a novel cooperative coupled particle filtering approach that enables us to track the dynamic changes of the shape, amplitude, and latency of the P300 subcomponents over different trials. This cooperative approach incorporates the spatial correlation of the ERPs and benefits from their bilateral symmetry to overcome the uncertainty of the estimation and tracking processes. Here, each one of the four relevant EEG channels acquire its own observation (the P300 waveform) to track the unknown underlying parameters. The channels can then share their information with each other in order to enhance the overall performance and obtain better estimates.

The proposed cooperative coupled Rao-Blackwellised particle filtering (CC-RBPF) is shown to be effective and robust for estimation of P300 subcomponents in single trials and can be applied to estimate other ERP subcomponents in single trials. This cooperative approach is novel and can be used in several applications dealing with simultaneous tracking of two dependent systems. The remainder of this paper is organized as follows. Section 2 introduces the particle filtering method and describes the Rao-Blackwellised particle filtering approach for estimating state variables. In Section 3, we first introduce the concept of coupled RBPF and then describe the proposed CC-RBPF method to estimate the state variables in a cooperative manner. In Section 4, the performance of the proposed method is evaluated and compared with the non-cooperative approach over several simulated data. Moreover, the method is applied to real EEG recordings from healthy and schizophrenic subjects in an oddball experiment. The paper is concluded in Section 5.

2. Particle Filtering

Particle filtering (PF) has been used as an effective approach to track the variations of P300 subcomponents over consecutive time-locked trials [32]. As discussed earlier, our goal is to track and estimate the values of amplitude, latency, and width of the subcomponents P3a and P3b among different trials. Therefore, we can form the state vector of the PF as:

$${\mathbf{x}}_k={[{\mathbf{a}}_k\left(1\right){\mathbf{b}}_k\left(1\right){\mathbf{s}}_k\left(1\right)\dots {\mathbf{a}}_k\left(p\right){\mathbf{b}}_k\left(p\right){\mathbf{s}}_k\left(p\right)]}^T$$

(1)

where ${\mathbf{a}}_k\left(i\right)$ represents the amplitude, ${\mathbf{b}}_k\left(i\right)$ is the latency, and ${\mathbf{s}}_k\left(i\right)$ is the width of the $i_{th}$ subcomponent in the $k_{th}$ trial, where $i=1,...,p$ and p is the number of subcomponents. Now, the state transition and observation equations can be formulated as:

$${\mathbf{x}}_k={\mathbf{x}}_{k-1}+{\mathbf{w}}_{k-1}$$

(2)

$${\mathbf{z}}_k=f(t;{\mathbf{x}}_k)+{\mathbf{v}}_k$$

(3)

where ${\mathbf{z}}_k$ represents the time-locked single trial measurement, ${\mathbf{w}}_{k-1}$ is the state noise, ${\mathbf{v}}_k$ is the measurement noise, and t is the time index which varies over the duration of the ERP component. Therefore, an ERP waveform can be modeled as:

$$f(t;{\mathbf{x}}_k)=f(t;{\mathbf{a}}_k,{\mathbf{b}}_k,{\mathbf{s}}_k)=\sum _{i=1}^p{a}_k\left(i\right)e^{-\frac{{(t-{\mathbf{b}}_k\left(i\right))}^2}{2{\mathbf{s}}_k^2\left(i\right)}},i=1,...,p.$$

(4)

It should be noted that the entire time samples are used in the observation resulting in a time vector [$f(1;{\mathbf{x}}_k),$ ... $,f($T$;{\mathbf{x}}_k)$], where T represents the number of time samples over the duration of the ERP component.

According to this model, it is assumed that each subcomponent can be approximated by a Gaussian waveform and the respective ERP component is formulated as the sum of these subcomponents. It has been shown that the subcomponents of P300 can be modeled with the help of parametric functions [33,34], among which the Gaussian waveform is the most common one [35]. Although the ERP subcomponents might not exactly follow the Gaussian or other functions, such modeling allows a fast and robust estimation to obtain the peak parameters of latency and amplitudes. The variations in these parameters are in fact the primary concern of cognitive scientists and neurophysiologists.

In order to track the ERP subcomponents, which consist of the amplitude, latency, and width of each subcomponent, and estimate the state of the system, the PF uses time-locked single trials as the available measurements. According to the PF approach, a posterior density function $p\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})$ is built to estimate the system state by a collection of random samples and their corresponding weights. This posterior density function is approximated as [36]:

$$p\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})\propto \sum _{n=1}^{N_s}{w}_k^n\delta ({\mathbf{x}}_k-{\mathbf{x}}_k^{\left(n\right)})$$

(5)

where $\{{\mathbf{x}}^{\left(n\right)},n=1,\dots ,N_s\}$ represents the particles with their associated weights $\{w^n,n=1,\dots ,N_s\}$. Since drawing samples from the unknown posterior density can be difficult, we employ the importance sampling approach to sample the particles from another density termed the importance density [37]. Therefore, the weights in (5) can be formulated as:

$$w_k^n=\frac{p\left({\mathbf{x}}_k^n\right|{\mathbf{z}}_{1:k})}{q\left({\mathbf{x}}_k^n\right|{\mathbf{z}}_{1:k})}$$

(6)

where $q(.)$ denotes the importance density function. In this paper, we employ the prior density function $p({\mathbf{x}}_k|{\mathbf{x}}_{k-1}^{\left(n\right)})$ as the importance density, which is the most common choice in the literature. Therefore, the weights can be updated by factorizing the importance density [36]:

$$w_k^n\propto w_{k-1}^{n}p({\mathbf{z}}_k|{\mathbf{x}}_k^{\left(n\right)})$$

(7)

It can be observed from (7) that the weight of each particle is updated according to its own weight in the previous trial and the observation likelihood by the corresponding particle given the current state. Over subsequent trials, the weights of a large portion of the particles might become quite small resulting in a very trivial contribution to the posterior estimation process while requiring significant computational cost for updating these weights. Therefore, we use the resampling technique which eliminates those particles that have small weights and replace them with particles with large weights. Moreover, initializing the particles effectively and placing them appropriately in the posterior density can be helpful to tackle this issue, especially for the cases with high dimensional state space.

2.1. Rao-Blackwellised Particle Filter

According to (4), the state space variables can be partitioned into linear and nonlinear variables, where the amplitudes are linear variables, and the latency and width comprise the nonlinear part. As a result, we can benefit from Rao-Blackwellised particle filters (RBPFs), where the state space is conditionally linear and some of the variables can be marginalized out analytically. In this way, the RBPF requires a smaller number of particles to provide the same performance as its counterpart PF.

In the RBPF approach, the state vector ${\mathbf{x}}_k$ can be partitioned as:

$${\mathbf{x}}_k=\left[\begin{array}{c}{\mathbf{x}}_k^1\\ {\mathbf{x}}_k^2\end{array}\right]$$

(8)

where ${\mathbf{x}}_k^1$ and ${\mathbf{x}}_k^2$ are the nonlinear and conditionally linear state variables respectively. To understand the process, we consider the case with the following linear Gaussian state space [38]:

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^1& \sim p\left({\mathbf{x}}_k^1\right|{\mathbf{x}}_{k-1}^1)\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^2& ={\mathbf{B}}_{k-1}\left({\mathbf{x}}_{k-1}^1\right){\mathbf{x}}_{k-1}^2+{\mathbf{w}}_{k-1}^2\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\mathbf{z}}_k& ={\mathbf{C}}_k\left({\mathbf{x}}_k^1\right){\mathbf{x}}_k^2+{\mathbf{v}}_k\hfill \end{array}$$

(11)

where ${\mathbf{x}}_k^1$ is a Markov process, ${\mathbf{B}}_{k-1}\left({\mathbf{x}}_{k-1}^1\right)$ is the state transition for the linear state variable, and ${\mathbf{C}}_k$ acts as a function of the nonlinear state variables. The function ${\mathbf{C}}_k$ can be combined with the linear state variable in order to model the observation. Moreover, ${\mathbf{w}}_{k-1}^2\sim N(0,Q_{k-1}^w)$ and ${\mathbf{v}}_k\sim N(0,Q_k^v)$ are zero mean Gaussian white noise distributions with known covariance matrices $Q_{k-1}^w$ and $Q_k^v$ respectively. In this paper, we select ${\mathbf{B}}_{k-1}$ as the unitary matrix, which makes the transition of the linear state variable independent of the nonlinear state variables. It can be observed from (10) that ${\mathbf{x}}_k^2$ has a linear Gaussian state space model subject to a given ${\mathbf{x}}_k^1$. We assume that the posterior density function $p({\mathbf{x}}_k^1,{\mathbf{x}}_k^2|{\mathbf{z}}_{1:k})$ can be reformulated as:

$$p({\mathbf{x}}_k^1,{\mathbf{x}}_k^2|{\mathbf{z}}_{1:k})=p\left({\mathbf{x}}_k^2\right|{\mathbf{x}}_k^1,{\mathbf{z}}_{1:k})p\left({\mathbf{x}}_k^1\right|{\mathbf{z}}_{1:k})$$

(12)

where we can analytically track and compute $p\left({\mathbf{x}}_k^2\right|{\mathbf{x}}_k^1,{\mathbf{z}}_{1:k})$ using Kalman filtering. Additionally, particle filtering can be utilized to estimate $p\left({\mathbf{x}}_k^1\right|{\mathbf{z}}_{1:k})$ [36]. Considering a conditionally linear Gaussian state space model, we can estimate the mean ${\mathbf{\Omega }}_k$ and covariance ${\mathbf{P}}_k$ of the linear state variable ${\mathbf{x}}_k^2$:

$$\begin{array}{c}{\mathbf{\Omega }}_k={\mathbf{\Omega }}_{k-1}+K_k({\mathbf{z}}_k-{\mathbf{C}}_k\left({\mathbf{x}}_k^1\right){\mathbf{\Omega }}_{k-1})\hfill \end{array}$$

(13)

$$\begin{array}{c}{\mathbf{P}}_k={\mathbf{P}}_{k|k-1}-{\mathbf{K}}_k{\mathbf{C}}_k\left({\mathbf{x}}_k^1\right){\mathbf{P}}_{k|k-1}\hfill \end{array}$$

(14)

$$\begin{array}{c}{\mathbf{K}}_k={\mathbf{P}}_{k|k-1}{\mathbf{C}}_k^T\left({\mathbf{x}}_k^1\right){({\mathbf{C}}_k\left({\mathbf{x}}_k^1\right){\mathbf{P}}_{k|k-1}{\mathbf{C}}_k^T\left({\mathbf{x}}_k^1\right)+{\mathbf{Q}}_k^v)}^{-1}\hfill \end{array}$$

(15)

$$\begin{array}{c}{\mathbf{P}}_{k|k-1}={\mathbf{P}}_{k-1}+{\mathbf{Q}}_{k-1}^w\hfill \end{array}$$

(16)

Now, we can tackle the problem of tracking P300 subcomponents using the RBPF. As shown in (4), we model each ERP component as the sum of its subcomponents, each represented by a Gaussian waveform. This observation model is then factorized into a linear part (amplitudes) and a nonlinear part (exponential waves). Using the RBPF concept here, the state vector is partitioned into linear and nonlinear parts:

$$\begin{array}{cc}\hfill & {\mathbf{x}}_k^1={[{\mathbf{b}}_k\left(1\right){\mathbf{s}}_k\left(1\right)\dots {\mathbf{b}}_k\left(p\right){\mathbf{s}}_k\left(p\right)]}^T\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & {\mathbf{x}}_k^2={[{\mathbf{a}}_k\left(1\right)\dots {\mathbf{a}}_k\left(p\right)]}^T\hfill \end{array}$$

(18)

$$\begin{array}{c}{\mathbf{x}}_k={\mathbf{x}}_{k-1}+{\mathbf{w}}_{k-1}\hfill \end{array}$$

(19)

$$\begin{array}{c}{\mathbf{z}}_k=\left[\begin{array}{c}{f}_1(t;{\mathbf{x}}_k^1)\dots f_p(t;{\mathbf{x}}_k^1)\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_k\left(1\right)\\ \dots \\ {\mathbf{a}}_k\left(p\right)\end{array}\right]+{\mathbf{v}}_k\hfill \end{array}$$

(20)

where

$$f_i(t;{\mathbf{x}}_k^1)=f_i(t;{\mathbf{b}}_k\left(i\right),{\mathbf{s}}_k\left(i\right))=e^{-\frac{{(t-{\mathbf{b}}_k\left(i\right))}^2}{2{\mathbf{s}}_k^2\left(i\right)}}fori=1,...,p$$

(21)

Based on the concept of RBPF, we partition the state vector into linear and nonlinear variables. As seen in (17) and (18), ${\mathbf{x}}_k^1$ is the nonlinear state vector containing the nonlinear variables (latency and width), and ${\mathbf{x}}_k^2$ is the linear part of the state vector containing the linear variables (amplitudes). Having these two linear and nonlinear parts, which are similar to those in (8) for RBPF, enables us to use RBPF for this problem. Moreover, the ERP waveform is formulated in (4), where the task (similar to the RBPF problem) is to estimate the values of $a_k\left(i\right)$, $b_k\left(i\right)$, and $s_k\left(i\right)$. Here, it can be observed that ${\mathbf{C}}_k$ in (11) is modeled as $f_i(t;{\mathbf{b}}_k\left(i\right),{\mathbf{s}}_k\left(i\right))$, which contains the nonlinear latency and width state variables.

3. Cooperative-Coupled Rao-Blackwellised Particle Filtering

Our objective is to better track the dynamic changes of the latency, amplitude, and width of the P300 subcomponents over different trials from multichannel EEG. To achieve this, we propose a novel cooperative-coupled Rao-Blackwellised particle filtering (CC-RBPF) method to incorporate the topographic information (locations) of P3a and P3b as well as their spatial correlation. In this cooperative method, each EEG channel is modeled as a node where each node aims to estimate and track the unknown parameters of the P300 subcomponents. In this section, we first formulate the tracking problem in a coupled manner and expand it to a cooperative PF approach which incorporates the bilateral symmetry of the ERPs.

3.1. Coupled Rao-Blackwellised Particle Filtering

As discussed in Section 2.1, the RBPF can be used to track the dynamics of amplitude, latency, and width of the ERP subcomponents over a single channel time-locked measurement. To obtain more accurate results, the averaged ERP of several trials is used to initialize the particles rather than generating random particles. However, noise and artifacts in the measurements can make it difficult to track the subcomponents accurately in the trials. This can result in the generation of invalid particles that do not satisfy some of the physiological constraints on the state variables. These constraints are in fact based on the existing knowledge of the features of the ERP subcomponents. By including them in our method, we can remove the invalid particles and eliminate their negative impact on the estimation process of the posterior density.

Some of these constraints can be directly applied to the estimates of a single channel. However, from a single channel it is not always possible to track the two temporally correlated subcomponents accurately, especially in the cases where there is a sudden change in one of the parameters of the waveforms. Therefore, applying more meaningful constraints to more channels and coupling the estimates can enhance the performance. As mentioned in Section 1, we track the dynamics of P300 subcomponents over four channels. The P3a and P3b subcomponents are considered to be spatially disjoint on the scalp while overlapping temporally [16]. One of the constraints on the subcomponents arises from the fact that P3a has a frontocentral distribution and is mainly toward the frontal region. Therefore, the amplitude of this subcomponent must be larger in the frontal channels compared to its counterpart in the parietal region. On the other hand, the amplitude of P3b is larger over the parietal channels compared to the frontal electrodes.

As seen in Figure 1, we have two symmetric electrodes in the frontal region (F3 and F4), and two symmetric electrodes (P3 and P4) over the parietal region of the brain. We formulate an RBPF for each of these channels. As the goal is to track the two subcomponents of P300 (P3a and P3b), the state variables for each RBPF can be simplified as:

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^1& ={\mathbf{R}}_k^i={\left[{\mathbf{b}}_k^i{\mathbf{s}}_k^i{\tilde{\mathbf{b}}}_k^i{\tilde{\mathbf{s}}}_k^i\right]}^T\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\mathbf{x}}_k^2& ={\mathbf{A}}_k^i={\left[{\mathbf{a}}_k^i{\tilde{\mathbf{a}}}_k^i\right]}^T\hfill \end{array}$$

(23)

where $i=1,2,3,4$ refers to RBPF1 to RBPF4 associated with channels F3, F4, P3, and P4 respectively. Moreover, ${\mathbf{b}}_k^i$ and ${\mathbf{s}}_k^i$ denote the latency and width of P3a over the $i_{\mathrm{th}}$ RBPF in the $k_{\mathrm{th}}$ trial while ${\tilde{\mathbf{b}}}_k^i$ and ${\tilde{\mathbf{s}}}_k^i$ represent the latency and width of P3b over the same RBPF and trial respectively. Similarly, ${\mathbf{a}}_k^i$ and ${\tilde{\mathbf{a}}}_k^i$ are the amplitudes of P3a and P3b respectively.

The state equations are then of the following form:

$$\begin{array}{c}\hfill {\mathbf{R}}_k^i={\mathbf{R}}_{k-1}^i+{\mathbf{w}}_{k-1}^{\rho i}\\ \hfill {\mathbf{A}}_k^i={\mathbf{A}}_{k-1}^i+{\mathbf{w}}_{k-1}^{Ai}\end{array}$$

(24)

where ${\mathbf{w}}_{k-1}^{\rho i}$ and ${\mathbf{w}}_{k-1}^{Ai}$ are zero mean Gaussian white noises with the covariance matrices ${\mathbf{Q}}_k^{w\left(\rho i\right)}$ and ${\mathbf{Q}}_k^{w\left(Ai\right)}$.

To couple the constraints over the RBPFs, we form two pairs, represented by ${\mathbb{N}}_j$ ($j=1,2$):

$${\mathbb{N}}_1=\{\mathrm{F}3,\mathrm{P}3\}$$

(25)

$${\mathbb{N}}_2=\{\mathrm{F}4,\mathrm{P}4\}$$

(26)

It is assumed that the subcomponents P3a and P3b are originally of the same shape among the channels in a pair ${\mathbb{N}}_j$ while having different amplitudes. Although this might not be exactly true in some real applications, coupling the two RBPFs is helpful in tracking the changes simultaneously where the signal specifications are uncertain due to inter trial variability. In this way, we can prevent sudden deviations from the true values in the estimation process. We can then assume that the widths of each subcomponent (P3a, and P3b) is the same over the pairs ${\mathbb{N}}_1$ and ${\mathbb{N}}_2$, coupling {RBPF1, RBPF3} and {RBPF2, RBPF4} with the same width. Regarding the latencies, a small difference is incorporated over the pairs for the estimation task. Hence, the same weight is assigned to the RBPF pairs for each subcomponent in each trial. On the other hand, the amplitudes are different and there is a small delay in the latencies of the subcomponents. Figure 2 shows the simulated noise free P300 waveform over the channels F4 and P4 in pair ${\mathbb{N}}_2$.

We can then generate the latencies and width of the subcomponents according to the prior density while estimating their amplitudes using Kalman filtering. As discussed earlier, some invalid particles might be generated during this process that do not satisfy the physiological constraints imposed on the method. Therefore, these particles are automatically removed and their weights are set to zero. The list of constraints imposed on the system are:

• The estimated amplitude of the P3a subcomponent at F (F3 and F4) channels should be larger than its counterpart in the P (P3 and P4) channels.
• The estimated amplitude of the P3b subcomponent at P (P3 and P4) channels should be larger than the amplitude of P3b in the F (F3 and F4) channels.
• As P3a has a relatively shorter latency than P3b, the estimated latency of P3a should be shorter than the estimated latency of P3b over the channels.

3.2. Cooperative Particle Filtering

As mentioned earlier, it is known that P300 has a symmetric distribution over left and right hemispheres of the scalp [28,29,30,31]. Incorporating this symmetric distribution and considering the spatial correlation of P300 in the system, motivates us to form a cooperative PF to enhance the accuracy of the tracking and estimation process. Cooperative approaches are of growing interest in problems with interconnected nodes due to their robustness, enhanced performance, and simplicity [39,40,41,42]. In these approaches, each node acquires its own observation (the P300 waveform here) and process it to track the unknown underlying parameters. Moreover, the connected nodes can share their information with each other in order to enhance their performance and obtain a better estimate.

Here, in our proposed CC-RBPF approach, each EEG channel is modeled as a node that tracks the unknown parameters (amplitude, latency, and width) of the P300 subcomponents. Due to the bilateral symmetry of the P300 subcomponents, the objective of the nodes on the right hemisphere (F4 and P4) is the same as their counterparts on the left hemisphere (F3 and P3 respectively). However, the signals of different nodes are distorted by independent noise processes that might be of different types or different strengths. Therefore, these nodes can cooperate with each other to suppress noise and obtain more accurate estimates of the unknown parameters.

According to Figure 3, neighboring nodes that are connected by a link can communicate with each other and their associated particles can share their latest estimates. Each particle can then linearly combine its own estimate with the received information from its neighbors to obtain a more accurate estimate. In other words, we have a collection of nodes (EEG channels) in this method that can cooperate and share information with each other. Therefore, two steps are followed:

• Each node (EEG channel) has its own observation (P300 waveform) and is interested to process it via the coupled-RBPF method described in Section 3.1 to obtain its unknown underlying parameters (amplitude, latency, and width). Moreover, we can apply some known existing constraints on the estimates obtained from the RBPF method to remove the invalid particles and their negative impact.
• Nodes can share their estimates with each other in order to obtain a better estimate and enhance their performance. The motivation behind this cooperation is the spatial correlation of P300, where it has a symmetric distribution over the left and right hemispheres of the brain. As a result, each particle in the RBPF method linearly combines its own estimate with the estimates of its neighbors in order to obtain a more accurate solution.

Before formulating the cooperation step among the agents it should be noted that the nodes can combine and allocate different weights to the information received from their neighbors. As shown in Figure 3, the weight that node i assigns to the information received from node ℓ can be denoted by $c_{\ell i}$. Moreover, the summation of the weights $c_{\ell i}$ over each node i should be equal to 1 and the weight node i allocates to its own estimate is $c_{ii}$. Therefore, the estimates are updated as follows:

$$\begin{array}{cc}\hfill {\left({\mathbf{R}}_k^n\right)}_i& =\sum _{\ell =1}^L{c}_{\ell i}{\left({\mathbf{R}}_k^n\right)}_{\ell }=\sum _{\ell =1}^L{c}_{\ell i}{\left[{\left({\mathbf{b}}_k^n\right)}_{\ell }{\left({\mathbf{s}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{b}}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{s}}}_k^n\right)}_{\ell }\right]}^T\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\left({\mathbf{A}}_k^n\right)}_i& =\sum _{\ell =1}^L{c}_{\ell i}{\left({\mathbf{A}}_k^n\right)}_{\ell }=\sum _{\ell =1}^L{c}_{\ell i}{\left[{\left({\mathbf{a}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{a}}}_k^n\right)}_{\ell }\right]}^T\hfill \end{array}$$

(28)

where $L=4$ is the number of channels, ${\left({\mathbf{R}}_k^n\right)}_i$ and ${\left({\mathbf{A}}_k^n\right)}_i$ are the estimates of particle n in channel i, and:

$$\sum _{\ell =1}^L{c}_{\ell i}=1.$$

(29)

It can be observed from (27) that the estimate of each particle n in each channel i of the nonlinear variables (latency and width) is the combination of its own estimate and all the nonlinear estimates ${\left({\mathbf{R}}_k^n\right)}_{\ell }$ received from the other channels. This means that particles are able to collaborate and share their information in order to obtain a more accurate estimate. This is also the same for estimating the linear variables in (28). Here, it is important to note that the weight $c_{\ell i}$ which is the weight node i assigns to the information received from node ℓ, should take into account the similarity of the two nodes. In other words, nodes with similar objectives and estimates should allocate higher weights to each other’s information while nodes with different objectives assign lower weight to that information. Therefore, the weights of the links connecting symmetric nodes {F3-F4} and {P3-P4} should be obviously higher than those of asymmetric pairs such as {F3-P4} and {F4-P3}.

After removing the invalid particles according to the existing constraints and updating the estimates, the weights of the remaining ones are normalized separately. In the coupled RBPFs, a particle pair is effective when both of its particles have relatively high weights. Therefore, in the cases where only the weight of one of the particles is high while the weight of its pair is relatively low, we decrease the weight for both of these particles. This is done by including a parameter in our method to update the weights of each particle pair together:

$$\begin{array}{c}\hfill {\left(w_k^n\right)}_{{\mathbb{N}}_j}^1={\left(w_k^n\right)}_{{\mathbb{N}}_j}^1\times {\alpha }^{\left(n\right)}\times {\lambda }^{\prime }\\ \hfill {\left(w_k^n\right)}_{{\mathbb{N}}_j}^2={\left(w_k^n\right)}_{{\mathbb{N}}_j}^2\times {\alpha }^{\left(n\right)}\times {\lambda }^{\prime }\end{array}$$

(30)

where ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^1$ is the weight of the particle in the first pair of ${\mathbb{N}}_j$, ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^2$ is the weight of its pair in ${\mathbb{N}}_j$, ${\lambda }^{\prime }$ is a constant parameter, and

$${\alpha }^{\left(n\right)}=1-|{\left(w_k^n\right)}_{{\mathbb{N}}_j}^1-{\left(w_k^n\right)}_{{\mathbb{N}}_j}^2|.$$

(31)

Then, in the resampling step, a scaled weight of the particle pairs is used as follows:

$${\left(w_k^n\right)}_{{\mathbb{N}}_j}=\lambda \times {\left(w_k^n\right)}_{{\mathbb{N}}_j}^1+(1-\lambda )\times {\left(w_k^n\right)}_{{\mathbb{N}}_j}^2$$

(32)

where $0\le \lambda \le 1$ is a constant parameter. The summary of the proposed CC-RBPF is given in Algorithm 1. In the next section, we evaluate the performance of the proposed CC-RBPF method and compare it with the performance of the non-cooperative approach.

Algorithm 1 Proposed Cooperative-Coupled RBPF Method

set $k=0$ and generate the coupled pairs ${\mathbb{N}}_1=\{F3,P3\}$ and ${\mathbb{N}}_2=\{F4,P4\}$

set $n=1,...,N_s$, where $N_s$ is the number of particles, and $i=1,...,4$ denotes RBPF1 to RBPF4 respectively.

Generate random numbers ${\left({\mathbf{R}}_0^n\right)}_{{\mathbb{N}}_j}^1$, ${\left({\mathbf{R}}_0^n\right)}_{{\mathbb{N}}_j}^2$ according to random uniform distribution considering the averaged ERP.

set ${({\mathbf{\Omega }}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^1$ and ${({\mathbf{P}}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^1=I$

set ${({\mathbf{\Omega }}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^2$ and ${({\mathbf{P}}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^2=I$

fork = 1 to $k_{max}$ ($k_{max}$is the number of trials, and $j=1,2$)

- generate random numbers ${({\mathbf{w}}_k^{\rho \left(n\right)})}_{{\mathbb{N}}_j}^1$ and ${({\mathbf{w}}_k^{\rho \left(n\right)})}_{{\mathbb{N}}_j}^2$

- set ${({\mathbf{R}}_k^n)}_{{\mathbb{N}}_j}^1={({\mathbf{R}}_{k-1}^n)}_{{\mathbb{N}}_j}^1+{({\mathbf{w}}_k^{\rho \left(n\right)})}_{{\mathbb{N}}_j}^1$ and ${({\mathbf{R}}_k^n)}_{{\mathbb{N}}_j}^2={({\mathbf{R}}_{k-1}^n)}_{{\mathbb{N}}_j}^2+{({\mathbf{w}}_k^{\rho \left(n\right)})}_{{\mathbb{N}}_j}^2$ satisfying the following constraints:

- ${\left(s_k^n\right)}_{{\mathbb{N}}_j}^1={\left(s_k^n\right)}_{{\mathbb{N}}_j}^2$ , ${\left({\tilde{s}}_k^n\right)}_{{\mathbb{N}}_j}^1={\left({\tilde{s}}_k^n\right)}_{{\mathbb{N}}_j}^2$

- $|{\left(b_k^n\right)}_{{\mathbb{N}}_j}^1-{\left(b_k^n\right)}_{{\mathbb{N}}_j}^2|\le \tau$ , $|{\left({\tilde{b}}_k^n\right)}_{{\mathbb{N}}_j}^1-{\left({\tilde{b}}_k^n\right)}_{{\mathbb{N}}_j}^2|\le \tilde{\tau }$

- update ${({\mathbf{\Omega }}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^1$, ${({\mathbf{\Omega }}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^2$, ${({\mathbf{P}}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^1$ and ${({\mathbf{P}}_k^{A\left(n\right)})}_{{\mathbb{N}}_j}^2$ for each particle using Equations (13)–(16)

- update the estimates according to the cooperative approach

- set the values of the weights $c_{\ell k}$

- ${\left({\mathbf{R}}_k^n\right)}_i={\sum }_{\ell =1}^L{c}_{\ell i}{\left({\mathbf{R}}_k^n\right)}_{\ell }={\sum }_{\ell =1}^L{c}_{\ell i}{\left[{\left({\mathbf{b}}_k^n\right)}_{\ell }{\left({\mathbf{s}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{b}}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{s}}}_k^n\right)}_{\ell }\right]}^T$

- ${\left({\mathbf{A}}_k^n\right)}_i={\sum }_{\ell =1}^L{c}_{\ell i}{\left({\mathbf{A}}_k^n\right)}_{\ell }={\sum }_{\ell =1}^L{c}_{\ell i}{\left[{\left({\mathbf{a}}_k^n\right)}_{\ell }{\left({\tilde{\mathbf{a}}}_k^n\right)}_{\ell }\right]}^T$

- calculate the weight of each particle using Equation (7)

- if ${({\mathbf{\Omega }}_k^{a\left(n\right)})}_{{\mathbb{N}}_j}^1\prec {({\mathbf{\Omega }}_k^{a\left(n\right)})}_{{\mathbb{N}}_j}^2$ or ${({\mathbf{\Omega }}_k^{\tilde{a}\left(n\right)})}_{{\mathbb{N}}_j}^1\succ {({\mathbf{\Omega }}_k^{\tilde{a}\left(n\right)})}_{{\mathbb{N}}_j}^2$ or ${\left(b_k^n\right)}_{{\mathbb{N}}_j}^1\succ {({\tilde{b}}_k^n)}_{{\mathbb{N}}_j}^1$ or ${\left(b_k^n\right)}_{{\mathbb{N}}_j}^2\succ {({\tilde{b}}_k^n)}_{{\mathbb{N}}_j}^2$

- set ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^1=0$

- set ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^2=0$

- end if

- normalize the weights ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^1={\left(w_k^n\right)}_{{\mathbb{N}}_j}^1/{\sum }_{n=1}^N{\left(w_k^n\right)}_{{\mathbb{N}}_j}^1$ and ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^2={\left(w_k^n\right)}_{{\mathbb{N}}_j}^2/{\sum }_{n=1}^N{\left(w_k^n\right)}_{{\mathbb{N}}_j}^2$

- set ${\alpha }^{\left(n\right)}=1-|{\left(w_k^n\right)}_{{\mathbb{N}}_j}^1-{\left(w_k^n\right)}_{{\mathbb{N}}_j}^2|$

- update the weights of both RBPF pairs

- set ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^1={\left(w_k^n\right)}_{{\mathbb{N}}_j}^1\times {\alpha }^{\left(n\right)}\times {\lambda }^{\prime }$

- set ${\left(w_k^n\right)}_{{\mathbb{N}}_j}^2={\left(w_k^n\right)}_{{\mathbb{N}}_j}^2\times {\alpha }^{\left(n\right)}\times {\lambda }^{\prime }$

- Resample particles of the pair RBPFs using the scaled weight

- set ${\left(w_k^n\right)}_{{\mathbb{N}}_j}=\lambda \times {\left(w_k^n\right)}_{{\mathbb{N}}_j}^1+(1-\lambda )\times {\left(w_k^n\right)}_{{\mathbb{N}}_j}^2$

end for

4. Experimental Results

In this section, we first apply the proposed method to simulated signals. The performance of this method is compared with the non-cooperative RBPF approach to evaluate its capability in estimating the P300 subcomponents. The method is then applied to real EEG data for healthy subjects and schizophrenic patients in an oddball paradigm.

4.1. Simulated Data

To evaluate the performance of the method, we generate two Gaussian waveforms representing P3a and P3b for each EEG channel. The sum of these two subcomponents at each channel forms the P300 ERP and the state variables to track are amplitudes, latencies, and widths of these Gaussian-shaped subcomponents. An example of the waveforms can be seen in Figure 2. To incorporate the bilateral symmetry of P300 in the simulated data, the noise free P300 waveform is the same over the symmetric channels {F3, F4} and {P3, P4}. The new state variables are generated from the previous trials using a uniform distribution, and the observations are obtained according to these new states. Moreover, the generated observations satisfy the physiological constraints on the latencies and amplitudes of the P3a and P3b subcomponents, and each subcomponent shares the same width with its counterpart in its coupled pair.

To more accurately model the real P300 ERP in the EEG, it is assumed that the P3a subcomponent has a shorter latency compared to P3b, and its amplitude in the frontal channels is larger than that of the parietal electrodes. In a similar manner, P3b has a larger amplitude over the parietal electrodes compared to the frontal ones. The simulations are run over 20 trials. Noise is added to the amplitude, latency and width of the subcomponents with variances 0.1, 3, and 1 respectively. Moreover, in each trial a random noise with variance 0.2 is added to the clean P300 observation, which is the summation of two Gaussian waves. The CC-RBPF method is then applied to the generated signals to track the P300 subcomponents over 20 trials. The values of the parameters $\lambda$ and ${\lambda }^{\prime }$ are empirically set as $\lambda =0.5$ and ${\lambda }^{\prime }=1$. Moreover, the values of $\tau$ (the difference between the P3a latencies on the two channels of a pair ${\mathbb{N}}_j$), and ${\tau }^{\prime }$ (the difference between the P3b latencies on the two channels of a pair ${\mathbb{N}}_j$) are set to 7. The number of particles used here is 10000. It is important to note that using RBPF rather than PF, allows us to decrease the number of required particles as it reduces the dimension of the state space. The reduced number of particles results in a lower computation time, especially in our method where several constraints should be checked for each particle.

The weights $c_{\ell i}$ in Equations (27)–(29) are designed in a way that the weight each node allocates to its own estimation is equal to the weight assigned to its symmetric node while the weight assigned to the other nodes is set as zero. For example, the weights assigned by node 1 (F4) in Figure 3 is set as: $c_{1,1}=c_{2,1}=0.5,c_{3,1}=c_{4,1}=0$. The sum of all the weight assigned by this node is equal to 1 and the weights of all the other nodes follow a similar pattern. As discussed earlier, it is crucial to design the weights such that nodes with similar observations allocate higher weight to each other. The performance of the proposed CC-RBPF method is also compared with a non-cooperative coupled RBPF algorithm, where the only difference is that in the non-cooperative case the symmetric nodes do not share their estimates and there is no cooperation among them in the estimation task. In other words, in the non-cooperative algorithm the weight that each node i assigns to its own estimate is $c_{i,i}=1$ while all the other weights $c_{\ell ,i}$ are set to zero.

Figure 4 and Figure 5 show the results of tracking amplitudes for the proposed CC-RBPF and non-cooperative approach respectively. In Figure 6 and Figure 7, the results of the obtained latencies are shown for CC-RBPF and non-cooperative approach respectively. Figure 8 and Figure 9 depict the tracking of the width values for CC-RBPF and non-cooperative methods respectively. As can be seen, the CC-RBPF approach can better estimate and track the variations of the P300 subcomponents. The existing deviations in the estimates are due to the random generation of the state parameters in each trial.

Figure 10 shows the signal to noise power ratio of each trial in dB defined as:

$$\mathbf{SNR}=10log\left(\frac{P_{signal}}{P_{noise}}\right)$$

(33)

Figure 11 presents the obtained mean square error (MSE) for both cooperative and non-cooperative methods, and is defined as:

$$\mathbf{MSE}=\frac{1}{N}\sum _{i=1}^N{(z_{act}\left(i\right)-z_{est}\left(i\right))}^2$$

(34)

where N denotes the number of time samples, $z_{act}$ is the simulated clean observation, and $z_{est}$ represents the estimated observation. The MSE obtained by the proposed CC-RBPF method is smaller then that of the non-cooperative approach for most of the trials in all the channels. These results show that the proposed cooperative method is able to track the inter-trial variabilities of the P3a and P3b subcomponents with a good accuracy. Its performance is also superior to the non-cooperative approach where the spatial correlation has not been taken into account.

4.2. Experimental EEG Data

The proposed method was applied to EEG recordings of healthy subjects and schizophrenic patients performing a two-stimulus auditory oddball experiment, previously used in [22]. The EEG data were recorded from 15 electrodes following the international 10–20 system with 2kHz sampling frequency. The subjects were asked to sit still with closed eyes to avoid any interference. Moreover, their necks were firmly supported by back of the chair to avoid muscle artifacts.

In the auditory oddball experiments, the stimuli are presented through earplugs, where an infrequent target occurs randomly in a frequent background stimuli. In this experiment, forty 1 kHz rare tones randomly occurred among 160 frequent 2 kHz tones. The intensity of these tones was 65 dB with duration of 10 ms and 50 ms for rare and frequent tones, respectively. The subjects were told to press a button when they heard the low target tone (1 kHz). All standard and target stimuli occurred every two seconds. The standard stimuli were consecutively repeated for 50 times before the first target occurs in order to increase the chance of P3a generation.

The proposed method was applied to the data from F3, F4, P3, and P4 channels considering the 40 rare trials, where both P3a and P3b are present. The average ERP of the first 20 trials from the total 200 trials was used to initialize the particles. The values of the parameters $\lambda$, ${\lambda }^{\prime }$, $\tau$, and ${\tau }^{\prime }$ are the same as those in the simulated data. Figure 12, Figure 13, Figure 14, Figure 15, Figure 16 and Figure 17 show the variability of P3a and P3b amplitudes, latencies, and widths over the rare trials for a healthy subject and a schizophrenic patient, respectively. It is observed from these figures that the amplitudes of P300 subcomponents of schizophrenic patients are much smaller compared to healthy subjects, which is in agreement with the existing literature [43,44,45,46]. In fact, one of the most consistent findings in schizophrenic patients is the significant reduction of P300 especially in auditory experiments [43,47]. Additionally, it can be observed that schizophrenic patients have longer P300 latency compared to healthy subjects, which is in agreement with the results of previous studies [22,43,44]. The mean width of P300 subcomponents is slightly larger for healthy subjects compared to schizophrenic patients. The method has been applied to eight more subjects and similar results have been obtained. According to these results, the amplitude of P300 subcomponents is a more discriminative feature for comparing schizophrenic patients with healthy subjects.

Based on the simulation results and experimental EEG dataset presented in this section, it is clear that the proposed CC-RBPF method is capable of detecting and tracking the dynamics of P3a and P3b subcomponents. The advantage of this method is that it can successfully estimate the P300 subcomponents in single trials and provides a more accurate and robust performance compared to the non-cooperative approaches. This is due to the cooperation among the nodes which helps them to suppress the existing noise in the system.

5. Conclusions

In this paper, we proposed a novel approach to detect and track the subcomponents of P300 over consecutive trials. We designed a cooperative coupled method using Rao-Blackwellised particle filters. With four electrodes, we formed two coupled RBPF pairs where a number of physiological constraints are imposed over each pair to obtain reliable estimations of the state variables. The neighboring channels share their estimates to obtain a more accurate estimation of the underlying parameters. This cooperative approach is motivated by the bilateral symmetry and spatial correlation of the ERPs over the scalp, helping the particles to suppress the noise and obtaining a more robust estimation. As a result, the proposed algorithm successfully estimates the P300 subcomponents in single trials and outperforms its non-cooperative counterpart. This cooperative approach can be effectively used for a variety of applications such as studying brain abnormalities or mental fatigue analysis. Also, the proposed approach can be utilized as a useful method to study the mental state of schizophrenic patients or other mental disorders. Since the method is capable of simultaneously tracking and characterizing two dependent systems, it can be also used for monitoring different natural systems. Another potential direction is the analysis of ERP dynamics over single trials using multivariate multiscale entropy analysis to gain a deeper insight into the complexity of ERPs and their variation over signal trials.