Separation and Localisation of P300 Sources and Their Subcomponents Using Constrained Blind Source Separation

Separation and localisation of P300 sources and their constituent subcomponents for both visual and audio stimulations is investigated in this paper. An e ﬀ ective constrained blind source separation (CBSS) algorithm is developed for this purpose. The algorithm is an extension of the Infomax BSS system for which a measure of distance between a carefully measured P300 and the estimated sources is used as a constraint. During separation, the proposed CBSS method attempts to extract the corresponding P300 signals. The locations of the corresponding sources are then estimated with some indeterminancy in the results. It can be seen that the locations of the sources change for a schizophrenic patient. The experimental results verify the statistical signiﬁcance of the method and its potential application in the diagnosis and monitoring of schizophrenia


INTRODUCTION
Event-related potentials (ERPs) are those electroencephalograms (EEGs) which directly measure the electrical response of the cortex to sensory, affective, and/or cognitive events. The fine-grained temporal resolution offered by ERPs allows accurate study of the time course of information processing unavailable to other neuroimaging techniques. However, spatial resolution has been traditionally limited. In addition, overlapping components of the ERP which represent specific stages of information processing are difficult to distinguish [1,2]. An example is the composite P300 wave, a positive ERP component which occurs with a latency of about 300 milliseconds after novel stimuli, or task relevant stimuli, requiring an effortful response on the part of the individual under test [1][2][3][4][5]. The P300 wave represents cognitive functions involved in orientation of attention, contextual updating, response modulation, and response resolution [1,3], and consists of multiple overlapping subcomponents, two of which are identified as P3a and P3b [2,5]. P3a reflects an automatic orientation of attention to novel or salient stimuli independent of task relevance [5,6]. Prefrontal, frontal, and anterior temporal brain regions play a major role in generating P3a giving it a frontocentral distribution [1,5]. In contrast, P3b has a greater centroparietal distribution due to its reliance on posterior temporal, parietal, and posterior cingulate mechanisms [1,2]. P3a is also characterised by a shorter latency and more rapid habituation than P3b [2,5]. Figure 1 illustrates some typical P3a and P3b waveforms from temporal-basal and temporo-superior dipoles [7].
Abnormalities in P300 are found in several of psychiatric and neurological conditions [4], however, differences may exist in particular in the specific subcomponents [2]. Moreover, changes to certain P300 subcomponents may distinguish between relatives discordant for psychiatric illness, and between subdiagnosis of illness [8,9]. That is, although reduced amplitude of the auditory P300 is reported in almost all studies of schizophrenia, the nature of these reductions including topography and associated subcomponents varies with subdiagnosis and sex [8,10]. Finally, certain subcomponents may be modality specific, whilst others may be independent of modality [2]. Thus, auditory and visual P300 appear to be differentially affected by illness and respond differently to treatment, suggesting differences in underlying structures and neurotransmitter systems [2]. P300 has significant diagnostic and prognostic potential especially when combined with clinical evaluation [2,4]. However, in order for this to be fully realised, efficient and reliable methods for separating P300 sources and its subcomponents must be established [4]. Blind source separation (BSS) has been used to identify the ERP subcomponents [11]. The objective of BSS is to separate a number of sources (component generators) from their mixtures (electrode signals). This is achieved by using information only from the sensor signals and, if available, some information about the statistical properties of the sources. Successfully performing BSS is a challenging problem in a variety of real-world applications. Various algorithms have been developed depending on specific applications [11]. A family of BSS algorithms stems from the principle of independent component analysis (ICA). This method tries to estimate the sources by assuming that they are statistically independent.
The most common method in detection, highlighting, and visualisation of P300 components used by clinicians is the frame averaging method. The problem has been tackled in more mathematical ways and one of the first approaches was to estimate brain sources, obtained from an electroencephalogram (EEG) or magnetoencephalogram (MEG), using a least-squares approach [12,13]. ICA was used later by a number of authors [14,15]. The motivation was to extract sources of electrical activity which represent different brain functions (i.e., they are independent). These authors used the Infomax [16] algorithm which produced satisfactory results in terms of source separation. In this paper, we develop a constrained algorithm based on Infomax to separate P300 sources and their subcomponents. The constraint term is achieved based on a prior knowledge of some measurable properties of the sources such as their latencies. A similar method has been developed in [17,18] which employs a constrained ICA algorithm using reference signals. However, the type of constraint and they way it is constructed are different. Here, we also emphasize on the development an automatic detection and localisation procedure for the P3a and P3b subcomponents.
The paper is structured as follows. Section 2 describes the basic BSS model and principles. Section 3 describes the proposed methods for separation and localisation of the P300 components. Section 4 shows some experimental results obtained by applying the proposed methods to EEGs recorded from some normal subjects and schizophrenic patients as well as some simulated data. Section 5 concludes the paper.

BLIND SOURCE SEPARATION
Regarding EEG, the mixing process is assumed instantaneous and the model is defined as follows: there is a number m of sensor signals x(t) = [x 1 (t), x 2 (t), . . . , x m (t)] T and a number n of source signals s(t) = [s 1 (t), s 2 (t), . . . , s n (t)] T . The mixing system is described by the matrix H ∈ R m×n and the relation between the sensors and the sources is The element h i j of the matrix H is the mixing coefficient from the jth source to the ith electrode. All signals are assumed zero mean or can be made so by subtracting the mean from them, additive observation noise is assumed insignificant. The aim of BSS is to estimate the original sources using information only from the sensors x(t). In most cases a matrix W = H −1 is used to estimate the sources indirectly by where s(t) denotes the estimate of s(t). In ICA the source signals are treated as random variables and the statistical properties of the signals are used to obtain the unmixing matrix. If each source i = 1, . . . , n is assumed to have a probability density function (pdf) q i (·), the independence assumption can be expressed mathematically as follows: the joint pdf q(s) of the source vector s is equal to the product of the marginal pdfs: The ICA algorithm usually depends on the assumption made for the pdfs of the sources q i (s i ). Also, since a closed-form solution to the ICA equation does not exist or it is generally very difficult to obtain, a cost function J(W), which provides a measure of independence, is optimised in an iterative manner using an optimisation technique such as a form of steepest descent or Newton's method. There are two main problems associated with the ICA method. Firstly, the estimated sources can be a scaled version (potentially with a sign change) of the original sources and, secondly, there is no way of knowing the order of the sources. These two problems are known, respectively, as the scaling and permutation ambiguities. The scaling problem may be mitigated by normalising the results with respect to the geometrical dimensions of the head. The permutation problem, however, has negligible effect in this application context as will be discussed in Section 4.

Constrained BSS
The Infomax algorithm was used as the original cost function since it has been reported to be effective for the separation Loukianos Spyrou et al.

3
of EEG signals [14,15]. The Infomax algorithm attempts to maximise the information flow between the inputs and the outputs of an artificial neural network (ANN). In this case, the inputs are the electrode signals and the outputs are some nonlinear transformation of the estimated sources. It is shown that if the nonlinear functions are selected appropriately [16], then the information maximisation will correspond to the minimisation of the dependence between the estimated sources. The Infomax cost function is where z ∈ R n×1 is the output of the neural network (z = f (y), f (·) is the nonlinear activation function applied element wise to y which is the estimated source vector), x is the input to the neural network, I(z, x) is the information between the inputs and the outputs of the ANN, H(z) is the entropy of the output and H(z | x) is the conditional entropy of the output assuming a known input; note, for convenience, the time index is dropped. The natural gradient of (4) is Maximisation based on the natural gradient is used to achieve good convergence [19]. The adaptation rule for the unmixing matrix W becomes where f (y)=(1 + exp(−y)) −1 with the assumption of super-Gaussian outputs and μ is the learning rate. The adaptation for an individual weight can be described by the equation (using the gradient ascent method) where cof represents the cofactor and det the determinant. Thus, each individual weight is adapted in a way that the rows and columns differ from each other, as prescribed by the first term of the right-hand side of the equation. When two rows or columns become similar, the matrix becomes singular, and then det W will tend to zero forcing the weight element to change dramatically. This change will be affected by cof w i j which shows the relative singularity of the remainder of the matrix, regardless of the row and column this element belongs to, compared to the whole of the matrix. ICA in general does not produce unique outputs and we aim to develop an algorithm that ensures that the desired P300 source is one of the estimated sources. This can be achieved by adding a constraint to the original algorithm. Lagrange multipliers incorporate the constraint function into the original cost function. This changes the problem into an unconstrained one. The constraint is considered as the Euclidean distance between the estimated sources and a reference P300 signal. The reference signal is obtained by frame averaging of the ERP obtained from a number of trials. The constrained problem can be written as J m and J C are the Infomax and the constrained cost functions, respectively. The cost function of the CBSS J algorithm is where Λ is the matrix of the Lagrange multipliers. The constraint function specialised for each column of W is defined as where r(t) is the reference signal and y i (t) is the ith output at time t. The unknown parameters in the problem are now two: the matrix W and the matrix Λ. The matrix W is found adaptively via the following relation [20]: where μ is the learning rate of the adaptation of the unmixing matrix, ρ is a scale factor for the Lagrange multiplier matrix, and P is a matrix whose rows contain the reference P300 signal. If a block algorithm is required, then the data vector x becomes a matrix and it should be scaled accordingly.
The basic form of the constrained algorithm can be modified to mitigate some inherent problems with this approach. Firstly, the present form of the algorithm tries to produce n outputs that are as close as possible to the P300 reference signal. Although this effect is alleviated partly by the Infomax algorithm which tries to produce different outputs, the constraint part of the algorithm will try and adjust more those outputs that are further away (in Euclidean distance terms) from the reference signal. Hence, it would be a good idea to try to enforce the constraint in one or a small number of the outputs. This comes from the fact that usually the P300 signal consists of a number of subcomponents in different regions of the brain. Secondly, the scaling ambiguity of every ICA algorithm can be a problem since one output could have exactly the same shape as the reference signal but it could be a scaled version of it. The algorithm would change that output (since it violates the constraint) which could damage its shape. So, a scaling procedure is used in which the reference signal matches the maximum amplitude of the estimated sources. Finally, the problem of finding good initial conditions for W, Λ, μ, and ρ, can be overcome partly by using a variable which determines the contribution of the two separate cost functions (i.e., main and constraint) to the adaptation of W. This way, the algorithm can be made to work (by avoiding the rapid divergence of the Frobenius norm of W) in a variety of situations. This way, the stability of the algorithm is ensured because the learning is kept bounded especially when J c (w i ) 0. It also functions as a safety point to make sure that the algorithm converges to a solution, which produces outputs close to the reference signals. The convergence of the algorithm is stable to the optimum point since both parts of the CBSS function have a negative definite Hessian matrix (easy to prove by checking the sign and the nonsingularity of the Hessian). The constrained cost function can take any form that would be suitable for a specific application. A cost function which maximises the inner product between the estimated sources and the reference signals was used but its performance was not as satisfactory as the Euclidean distance function. Following the theory of constrained optimisation, in cases where the separation needs to be improved over the traditional ICA methods, a number of new BSS algorithms can be developed. Other suggested cost functions for the present purpose can be maximising the spikiness of the output sources around the time of interest (300 milliseconds), estimating the pdf of the P300 sources and forcing the pdfs of the output sources to have a similar form 1 or even applying a spatial constraint using prior knowledge of the possible P300 positions.
A variation of this algorithm which was used to separate the P3a and P3b subcomponents was implemented by using the method of least squares. If the reference signals for P3a and P3b are known, then we can model the EEG system as where r is the reference signal, X is the data matrix, and w opt (row vector 1 × n) is the vector that should produce r. Then, the constraint cost function will be where w i is the ith row vector of the unmixing matrix. This vector corresponds to the ith output y i expected to be the separated P3a or P3b. The selection of the appropriate y i to enforce the constraint is achieved in terms of which one is closer in terms of the Euclidean distance to the reference signal. w opt is found using the common least-squares (LS) solution: Also, the gradient of (14) is Then, this gradient is incorporated within the main Infomax update equation in a similar manner to (11). This constraint is different from those used in [17,18].

Construction of the reference signals and detection of P300 subcomponents
P3a and P3b are the two P300 subcomponents that overlap at the scalp. A constrained BSS algorithm such as that described above can be used to extract the P3a and P3b from multichannel EEGs. One important factor in applying CBSS is the selection of the proper reference signal. The way we obtain the reference signals is to use prior knowledge of the latencies of the two subcomponents. P3a peaks on average at a latency of 260 milliseconds and P3b on average at 300 milliseconds. However, it is possible that both the P3a and P3b occur with different latencies. The distinctive feature is then that P3a occurs before the P3b. P3a is hence selected by space-time averaging all the electrodes and selecting the first peak that occurs near the time of interest (250 milliseconds-350 milliseconds) and P3b by selecting the second peak. The two reference signals are then used in the CBSS algorithm. To detect which of the CBSS outputs is the P3a and which is the P3b, we use the correlation function. For two variables x and y, the correlation coefficient is defined as where σ x denotes the standard deviation and cov(x, y) the covariance of the two variables. The covariance of the two variables provides a measure of how strongly correlated these variables are. Because our purpose is to detect P3a and P3b, the source which has the maximum correlation coefficient with the P3a or P3b reference signal is more likely to be P3a or P3b, which will be selected automatically.

Localisation
Localisation of electrical sources inside the brain has been investigated by a number of people [12,13,[21][22][23]. Unlike the work done in [12] which assumes that the electrical sources are magnetic dipoles, here we assume that they are sources of isotropic propagation. Hence, the head simply mixes and attenuates the signals. Therefore based on Figure 2 we have where f k is the position of the source k, a j are the positions of the electrodes, and d j are the distances between the source and the jth electrode. The distances d j are nonlinearly proportional to the inverse of the correlation between the estimated source and the electrode signals. This is because a source is attenuated nonlinearly with the distance. Hence, the correlation of the electrodes with a source is nonlinearly proportional to the distance [24]: where X = HS, H describes the forward model for which the magnitude of a source attenuates with 1/d 2 j , and s k is the vector of all sample values of the kth source. It has to be noted that the sources must be uncorrelated for the method to be efficient. After computing the correlation the values are normalised and converted to distances by the following: It has to be noted that this approach does not provide a valid source reconstruction since it ignores the conductivity properties of the brain but it can be used to distinguish between sources in relatively different locations. Index Loukianos Spyrou et al. The next step is to convert to a mathematical problem, which is required to calculate the coordinates of an unknown point when both of the coordinates of d points and the distances of the unknown point from the given points are known. This problem is clearly equivalent to finding the intersection point(s) of d spheres in R d . The points are the solution to the following least-squares problem which can be obtained from [25]: where

Simulated data experiment
The CBSS algorithm was firstly applied to simulated data to test its efficacy. Two sinusoidal sources and one sinc signal (simulating the P300 as in Figure 3) were mixed (Figure 4). The results from the original Infomax algorithm were compared with those of the CBSS algorithm. The same learning rate (μ = 0.001) and initial conditions (W 0 = I) were used in both cases. Figures 5 and 6 show, respectively, the estimated P300 source using the two algorithms. The original sinc source was used as the reference. The performance was measured with the mean square error: where N is the number of samples, s(i) is the obtained source, and s(i) is the original source and the error is given in dB. The P300 obtained from normal Infomax had an error of e = −18.9 dB while that of CBSS had an error of e = −19.8 dB.
The CBSS algorithm used was that of (11).

Experiment data
The EEG data were recorded using a Nihon Kohden model EEG-F/G amplifier and Neuroscan Acquire 4.0 software. EEG activity was recorded following the international 10-20 system from 15 electrodes. The reference electrodes were linked to the earlobes. The impedance for all the electrodes was below 5 kΩ, sampling frequency Fs = 2 kHz, and the  Figure 6: The estimated P300 sources using Infomax. The simulated P300 is slightly more distorted than that obtained by CBSS and achieves a higher mean square error.
data were subsequently bandpass filtered (0.1-70 Hz). This frequency range was chosen to be compatible with [26]. The EEG data were recorded for control and schizophrenic patients. The Diagnostic and Statistical Manual of Mental Disorders 4th edition (DSM-IV) Axis 1 disorders was used to confirm diagnosis of schizophrenia. Subjects were required to sit alert and still with their eyes closed to avoid any interference. Also, to avoid any muscle artefact, the neck was firmly supported by the back of the chair. The feet were rested on a footstep. The stimuli were presented through ear plugs inserted in the ear. Forty rare tones (1 kHz) were randomly distributed amongst 160 frequent tones (2 kHz). Their intensity was 65 dB with 10-and 50-milliseconds duration for rare and frequent tones, respectively. The subject was asked to press a button as soon as they heard a low tone (1 kHz). The ability to distinguish between low and high tones was confirmed before the start of the experiment. The task is designed to assess basic memory processes. ERP components measured in this task included N100, P200, N200, and P3a and P3b.

Separation of P3a and P3b
Firstly, the ERP is obtained by temporally averaging eventrelated data (40 events), each event producing an EEG of size n × T, where n is the number of electrode signals and T is the number of samples of the event. That averaged ERP is also of dimensions n × T. The advantage of averaging event-related data is not only to enhance the signal, but also to remove non-event-related noise. Secondly, the reference subcomponent signal is selected according to the method described in Section 3.2. Thirdly, CBSS is applied to the ERP (n × T) in order to separate the P300 and its sub-components. Filtering (at the Delta range) is applied to the separated sources, based on the knowledge that the main power of the P300 component is in the Delta range [27]. Figure 7 shows the ERP and Figure 8 shows the estimated P3a and P3b sources for a schizophrenic patient. Figure 9 shows the ERP and Figure 10 shows the estimated P3a and P3b for a control subject. It can be seen that the P3a component is earlier in latency than the P3b.

Localisation of P3a and P3b
To approximately specify the location of a source in the head, we consider a spherical model of the head and assume isotropic propagation of the sources. Using the method Loukianos Spyrou et al. described in Section 3, there may be some trivial solutions (i.e., the points which fall outside the head) which are automatically discarded based on geometrical constraints. The result of the localisation of the P3a and P3b components is shown in Figure 11 for five schizophrenic patients and in Figure 12 for five control subjects. It is evident that the P3a and P3b for a schizophrenic patient are closely and irregularly located, whereas for a control subject the P3a and P3b are located in distinct regions.

Comparison between CBSS and Infomax
Some obtained P3a using normal unconstrained Infomax and CBSS is shown here. The results from normal Infomax, CBSS, and the reference signal are shown in Figures 13 and  14. It is seen that the CBSS produces better results than those of the normal Infomax in terms of highlighting of the relevant signals (P3a's latency is about 260-280 milliseconds). In quantitative terms CBSS can produce results with up to 33% more similarity with the reference signal. 2 Another significance of the proposed CBSS algorithm is that the P3a and P3b are robustly extracted. While Infomax may fail to produce those outputs (due to the nonstationarity of the data and the initialisation procedure), CBSS ensures that the desired outputs are always extracted.

Visual and auditory P300 comparison
The approximate source localisation method described in Section 3 was implemented for audio and visual ERPs separately. This was done to examine any differences in the locations to be further used in diagnosis of the psychiatric disorders. A set of EEG was obtained using the same hardware and software but with a 64-electrode cap.
To obtain the visual P300 the experiment consisted of a series of letters displayed successively with a period of 5 seconds. The image lasted 100 milliseconds. When a letter was displayed twice in a row, the subject had to press a button,  which should elicit a P3b. Occasionally, a checkerboard was displayed on the screen resulting in a P3a. The experiment lasted about 7 minutes. A similar experiment was performed to obtain the auditory P300. The sounds of different letters were played through ear plugs inserted into the ear. A sequence of letters was pronounced and when two were pronounced in series, the subject had to press a button, which should elicit a P3b. Intermittently, noise sound was played resulting in a P3a. The period was again 5 seconds. The experiment lasted about 7 minutes. The data which should elicit a P3b were selected for this experiment. CBSS was used to extract the P300 based on (11) and then the inner-product between the estimated P300 source and each electrode was computed. The estimated P300 from visual stimuli is shown in Figure 15 and from auditory stimuli in Figure 16. Figures 17 and 18 show the inner-product of the estimated P300 and the electrode for visual and auditory Loukianos Spyrou et al. Amplitude (mV) Figure 16: Auditory P300 obtained using CBSS. Figure 17: Distribution of visual P300 over the scalp electrodes. It can be seen that the P300 is distributed in a different way over the electrodes than the auditory P300. Figure 18: Distribution of auditory P300 over the scalp electrodes. It is seen that the auditory P300 is distributed differently over the electrodes than the visual P300. data, respectively. Results from two data frames are shown. It is observed that the latency for the visual P300 is longer than for the auditory P300. Another important conclusion is that the projections of the P300 audio and visual sources over the electrode are different. This means that these components are generated in different regions of the brain.

CONCLUSIONS
In this paper, a constrained BSS method has been developed to separate and localise the P300 signals and their constituent subcomponents from the EEG/ERP signals. The incorporated constraint minimises the distance between a measured reference signal and the estimated independent components. The proposed CBSS method achieves better performance in terms of extraction of the relevant signals. The algorithm was applied for separation and localisation of both audio and visual P300 sources. The CBSS method was also used to separate and localise the P3a and P3b subcomponents. A number of experiments on healthy subjects and patients suffering from schizophrenia were carried out. As a result, the latency of the P300 for schizophrenic patients was seen to be longer than that of the healthy person. Also, it was concluded that the P3a and P3b subcomponents are often located in completely different regions of the brain for the healthy subject whereas for the schizophrenic patients the sources are closely and irregularly located. Although the localisation algorithm has yet to be modified to mitigate the indeterminacy and to incorporate the nonhomogeneity of the head, the primary outcomes of this work are very valuable for diagnosis, treatment, and monitoring certain psychiatric illnesses.