Neural Computation, 27(7):1373–1404, 2015. This is a preprint version. The final publication
is available at http://www.mitpressjournals.org/doi/abs/10.1162/NECO_a_00747

Unifying Blind Separation and Clustering for Resting-State
EEG/MEG Functional Connectivity Analysis
Jun-ichiro Hirayama1,* , Takeshi Ogawa1 , and Aapo Hyvärinen2,1
1

Brain Information Communication Research Laboratory Group, Advanced
Telecommunications Research Institute International (ATR), Kyoto, Japan.
2
Department of Computer Science and HIIT, University of Helsinki, Helsinki, Finland.
*hirayama@atr.jp

Abstract
Unsupervised analysis of the dynamics (non-stationarity) of functional brain connectivity during rest has recently received a lot of attention in both the neuroimaging and neuroengineering communities. Most studies have used functional magnetic resonance imaging(fMRI), but electroencephalography (EEG) and magnetoencephalography (MEG) also holdgreat promise for analyzing non-stationary functional connectivity with high temporal resolution. However, previous EEG/MEG analyses divided the problem into two consecutive
stages: first, the separation of neural sources, and second, the connectivity analysis of theseparated sources. Such non-optimal division into two stages may bias the result because of
the different prior assumptions made about the data in the two stages. Here, we propose a
unified method for separating EEG/MEG sources and learning their functional connectivity(coactivation) patterns. We combine blind source separation (BSS) with unsupervised clustering of the activity levels of the sources in a single probabilistic model. A BSS is performedon the Hilbert transforms of band-limited EEG/MEG signals, and coactivation patterns arelearned by a mixture model of source envelopes. Simulation studies show that the unified
approach often outperforms conventional two-stage methods, further indicating the benefit
of using Hilbert transforms to deal with oscillatory sources. Experiments on resting-state
EEG data, acquired in conjunction with a cued motor imagery/non-imagery task, also showthat the states (clusters) obtained by the proposed method often correlate better with physiologically meaningful quantities than those obtained by a two-stage method.

1

Introduction

Unsupervised machine learning techniques play a fundamental role in the analysis of spontaneous(resting-state) neuroimaging signals by exploring the intrinsic statistical structures of such signals without relying on extrinsic covariates about tasks or stimulation protocols. The structures
or features obtained can then be examined based on neurophysiological knowledge often using
the features in a group comparison, or possibly by finding similar structures in other signals
already associated with tasks or stimuli.
In recent years, there has been growing interest in exploring the patterns and dynamics
of resting-state functional brain connectivity [11] based on unsupervised signal analyses. To
find patterns in non-stationary functional connectivity, most studies have relied on such standard techniques as independent component analysis (ICA) [6, 45], principal component analysis(PCA) [31], and K-means clustering [32, 1]. New unsupervised analysis methods have also1been actively developed [e.g., 16, 50, 25, 17, 39, 40, 8] to incorporate the specific nature ofneuroimaging signals.
Such methods for analyzing functional brain connectivity, primarily based on signals’ own
statistics, are potentially very useful not only for neuroscientific investigations but also forneuroengineering applications. For example, a key challenge in an emerging new direction inbrain-computer interface (BCI) research is to covertly acquire user’s unobserved states during
everyday life behaviors [49, 30]. Since brain activity cannot be very well controlled in everyday
life situations, no reliable class labels are available for discriminating the unobservable statesof users. Such difficult data could be tackled if unsupervised analysis discovered connectivity
patterns reflecting the user’s cognitive states. Unsupervised connectivity analysis might also
shed light on the neurophysiological basis of the BCI paradigms commonly used so far, such as
those based on motor imagery, i.e., imagining body movements [14, 15].
Motivated by such potential applications, in this paper we focus on developing an unsupervised analysis method for finding connectivity-related signal features from electroencephalography (EEG); our method may also be readily applied to magnetoencephalography (MEG) becauseof its fundamental similarity. EEG/MEG’s high temporal resolution is particularly useful for
analyzing non-stationary connectivity, as compared to functional magnetic resonance imaging(fMRI), which is used in most connectivity studies.
The analysis of functional connectivity patterns in EEG/MEG, however, is not straightforward because the neural sources are mixed by volume conduction (and/or field spread) into
sensor signals. Two-stage analysis has been conventionally performed by first separating the neural sources from the given sensor signals and analyzing the connectivity patterns based on thoseseparated sources. In neuroimaging literature, electromagnetic inverse problems are often solved
to separate (estimate) the activity of dipolar sources on the cortical grid, with an additional
effort of physical forward modeling. On the other hand, in exploratory signal analysis related
to BCI, blind source separation (BSS) methods (including ICA) are especially useful, since theygreatly simplify the interpretation of the results by decomposing the data into components, andthe inverse problem can be solved to localize the obtained components afterwards [19, 10]. Notethat these components are actually called “sources” in BSS literature, and a component can be
an integration of multiple correlated dipolar sources.
The problem is that conventional two-stage analysis, i.e., first source separation and thenconnectivity analysis, is “neither a principled nor an optimal solution to the overall problem” [34].
Source separation methods rely on specific prior assumptions about the sources, which are notnecessarily consistent with what the connectivity analysis assumes about them; such inconsistency between prior assumptions might bias the results. A more desirable unified treatment
would be obtained by extending the conventional generative (forward) model used for source
separation by including prior assumptions about the sources that are consistent with the connectivity analysis. Typically, connectivity analysis can be formulated as learning of a specificparametric model of the sources, and both layers of the model may be learned simultaneously,
unified by the principle of statistical parameter estimation.
Here, we present a unified method for BSS and analyzing functional connectivity patterns
in EEG/MEG sources with a jointly solved BSS, based on a novel two-layer extension of theconventional BSS/ICA generative model. In line with previous resting-state MEG studies [9, 6,40], we are particularly interested in finding coherent (frequently occurring) patterns of activity
levels (coactivations) of oscillatory sources in any frequency band of interest. The connectivities
are based on envelopes and ignore phase information, but our model is rather different frompower-to-power coherences.
To properly model the source envelopes, our model uses a complex-valued formulation of
BSS based on the Hilbert transform, which is a key departure from related extensions of ICA

2based on modeling real-valued data [21, 22, 47, 26, 27, 37, 29, 16, 17]. Another importantnovelty here is using a finite mixture model of sources, in which they are assumed to exhibit
different coactivation patterns corresponding to a finite number of unobserved “brain states.”
This corresponds to performing unsupervised clustering on the coactivations of the sources,
inspired by the use of K-means clustering in previous resting-state fMRI studies [32, 1]. Due tothe simplicity of the mixture model, our two-layer model is tractable, unlike previous two-layer
models that are often difficult to learn. Our entire model can be readily estimated (optimized)by the maximum likelihood method without resorting to any approximations.
The rest of this paper is organized as follows. First, we present our proposed method basedon a novel two-layer extension of the generative BSS/ICA model called the latent coactivity
mixture model (LCMM) (Section 2). Then we provide simulation studies (Section 3) and real
EEG data analysis (Section 4) to validate the unified approach. Finally, we discuss the results
and open issues (Section 5). Preliminary results were presented at the International Conferenceof the IEEE Engineering in Medicine and Biology Society (EMBC’14) [18].

2

Latent coactivity mixture model

2.1

Background: blind source separation

Before introducing our new model, we start by discussing the generative (forward) model conventionally used for the blind separation of EEG/MEG sources. Let x(t) = (x1 (t), x2 (t), . . . , xd (t))T ∈Rd be a multivariate EEG/MEG signal, sampled at discrete time points indexed by t = 1, 2, . . . , N ,which is assumed to have already been (band-pass) filtered so that each xj (t) is limited to acertain frequency band of interest. Sensor signal vector x(t) is then assumed to follow a linear
generative model given byx(t) = As(t),(1)

where s(t) = (s1 (t), s2 (t), . . . , sd (t))T ∈ Rd denotes the vector of the source signals and A ∈ Rd×d
is called the mixing matrix, which is assumed to be non-singular so that demixing matrix
W := A−1 exists. Both s(t) and A are unknown and estimated from data in the BSS setting.
Note that in Eq. (1), the number of sources is assumed to equal (effective) dimensionality
d of the sensor signal, as a fundamental setting in a standard ICA; this greatly simplifies themathematical treatment. In practice, d can be selected to be smaller than the original numberof sensor channels, typically by discarding the ineffective dimensions (with too small variances)using PCA. It should also be noted that each source (or component) sj (t) does not necessarily
correspond to any single electrical dipole; instead, each sj (t) may describe the total effect ofmultiple correlated dipolar activities.
To solve the BSS problem, we need to make further assumptions about the statistical properties of the sources based on prior knowledge. Independent component analysis (ICA) typicallyassumes that d sources are non-Gaussian and mutually independent [23], which theoreticallyguarantees the identifiability of both A and s(t)’s, up to the scaling and permutation of thesources. However, the independence assumption might lead to a solution that is weakly functionally connected (i.e., weakly statistically dependent), even when the true sources are strongly
dependent, which is not consistent with the goal of connectivity analysis. This motivated us todevelop an appropriate model of functionally connected (dependent) sources.

2.2

Definition of latent coactivity mixture model

Our main interest here is modeling the envelopes (i.e., amplitudes) of narrow-band source signals[36, 50, 6] for which the “real-valued” BSS model of Eq. (1) is not convenient. Envelopes can
3be modeled more easily with complex analytic signals sej (t) [43], defined assej (t) = sj (t) + iH[sj ](t),(2)

where H[f ] denotes the Hilbert transform of signal f (t) and i is an imaginary unit. The envelopeof sj (t) is given by the modulus of sej (t). This simple algebraic dependence of the envelope on se
greatly simplifies the developments below.
We thus formulate a “complex-valued” BSS with a similar transformation of sensor signal
vector x(t):e (t) = As
e(t),
x

(3)

e (t) = (xe1 (t), xe2 (t), . . . , xed (t))T ∈ Cd can be directly computed
where complex sensor signal x
ej (t) := xj (t) + iH[xj ](t). Here, sources sej (t) are further assumed to befrom the original one by x
centered (i.e., E[sej (t)] = 0) without loss of generality by always subtracting the (sample) meane . Both mixing matrix A and complex-valued source signal se(t) = (se1 (t), se2 (t), . . . , sed (t))Tfrom xare again unobserved and estimated from the data.
Note that mixing matrix A in Eq. (3) is identical to the original one in Eq. (1) because of
the Hilbert transform’s linearity. This ensures that the A’s columns can be interpreted directly
as defining spatial topographies in the original sensor space. For simplicity, we constrain A to
be real-valued, while a complex-valued A could also be straightforwardly used, which might beuseful to deal with sources synchronized in different phases [24].
We next define a coactivation (connectivity) structure between sources sej (t). The fundamental assumption here is that a system generating the data can be in a finite number of different
states, corresponding to different patterns of source amplitudes. Given the state at time pointt, sources sej (t) are generated based on a multivariate Student-t distribution (specified below)
which implements the average source amplitudes specific to that state and generates random
phases.
Thus, our model of connectivities is not based on explicitly measuring some form of correlations between the sources or their envelopes, as is typically done in electrophysiology. Instead,
we characterize the interactions of the sources by dividing their joint activity into a numberof typical patterns of envelopes, which intuitively express the idea that certain sources tend
to be coactivated. Such coactivation does imply correlations of envelopes, i.e., power-to-power
coherence, but provides a more detailed analysis of the coactivation than merely computing
correlations.
The proposed latent coactivity mixture model (LCMM) is thus summarized as a two-layere (t) as follows:
generative model of complex sensor signal vector x

1. At each time point t, the system generating the data takes one of a finite number ofdifferent states (clusters) indexed by k = 1, 2, . . . , K, according to multinomial probability
P
distribution with cluster probabilities η1 , η2 , . . . , ηK , where ηk ≥ 0 and K
k=1 ηk = 1.
e(t) is specifically
2. Given that the system belongs to the k-th state at time t, source vector s
generated by a complex multivariate Student-t distribution [43] with circular (see below)
and mutually uncorrelated sources; they have state-conditional variances or expected powers (squared amplitudes), given by

E[|sej (t)|2 ]k =

ν
bjk ,
ν−2

j = 1, 2, . . . , d,(4)where E[·]k denotes the conditional expectation given the k-th state, ν is an integer called

4the degrees of freedom, and nonnegative vector bk = (b1k , b2k , . . . , bdk )T , called the coactivation pattern, specifies the expected levels of the source envelopes.1
e (t) is given as a linear instantaneous mixture of se(t) with unknown
3. Complex sensor signal x
mixing matrix A, as in Eq. (3), common to each state k.
e (1), xe (2), . . . , xe (N ), are simply assumed to
The N Hilbert-transformed sensor signal vectors, x
be independently and identically distributed (i.i.d.), as is commonly done in many ICA methods.
Such a simplification is mainly done for purposes of mathematical and computational tractability,
but it could be relaxed by further modeling the autocorrelation structures of the sources (infuture work).
The circularity of the sources means that the phase of each source is distributed uniformly
and independently of its amplitude. Our method focuses on modeling (analyzing) an amplitudeto-amplitude type of connectivity (coactivation), ignoring phase-to-phase or phase-to-amplitude
types of connectivity. In the circular case, the probability density function (pdf) of complex
multivariate Student-t distribution [43] with scatter matrix Σ and ν (> 0) degrees of freedom
is given byν

2d Γ(d + ν2 )2 H −1 −d− 2
e Σ s
e
,
1
+
s
Γ( ν2 )(νπ)d |Σ|νe; Σ, ν) =
Te (s



(5)

0
~]
Re[s
1

1

(a) e
s1 : Real and imaginary parts

2

−2

(b) e
s1 and es2 : given a state

2

0

1

2
−1
0
~ ] or
~1 ]
Re[s
Im[s
1
1

−1

−1

0

1

~ ] or Im[s
~]
Re[s
2
2

2
−1

−2

−2

−2

0
−2

−1

~]
Im[s
1

1

~ ] or Im[s
~]
Re[s
2
2

2

where ·H denotes the Hermitian transpose and Γ(·) denotes the Gamma function. Source vector
e in LCMM has a state-conditional pdf given by Te (se; diag(bk ), ν). Figs. 1 (a) and (b) illustrates
the conditional pdf for a bivariate case.

−2

−1
0
~ ] or
~1 ]
Re[s
Im[s
1
1

2

(c) Mixture density

Figure 1: Illustration of complex multivariate Student-t distribution for two circular complex
sources. (a,b): Bivariate Student-t pdf Te ((se1 , se2 ), diag(b1 ), ν) with b1 = (1, .1) and ν = 2;(a) illustrates pairwise marginal density on real and imaginary parts of single source se1 , where
spherical equiprobability contours imply circularity; pairwise marginals for the four combinationsof real or imaginary parts between se1 and se2 have the same form illustrated in (b) (blacksolid lines), where two other examples are also shown (gray solid lines) with b2 = (.1, 1) andb3 = (.8, .8). (c): Mixture density of the three state-conditional pdfs in (b), with the state kmarginalized out with η = (.1, .1, .8).
1

To be precise, the exact relation between bk and the source variances in Eq. (4) is no longer valid if ν ≤ 2,
since the variance is infinite or undefined. However, even if ν ≤ 2, the estimated bk can still be interpreted asmodeling the variance levels, while small ν implies canceling the effects from the outliers from the viewpoint ofrobust estimation (as explained in the text).

5

(a) State k marginalized(b) With finite ν(c) With ν → ∞ (Gaussian)Figure 2: Graphical representations of dependency structure in LCMM
The particular choice of the Student-t model is mainly motivated by its robustness to outliersin the estimation of (co)variances (i.e., coactivation patterns), as has been thoroughly studied
in the literature [see e.g., 35, 33].2 Scatter matrix Σ is proportional to the covariance matrix ifeseH ] = {ν/(ν − 2)}Σ for ν > 2, and maximum likelihood estimate (MLE) of Σ
it exists, i.e., E[sis often used as a robust alternative of the sample covariance matrix even with ν ≤ 2. Diagonal
scatter matrix diag(bk ) in LCMM implies Eq. (4), and the estimate of bk serves as a robust
estimator of the state-conditional variances; note that even if the variance does not exist whenν ≤ 2, bk can have a finite MLE, giving a more general scale parameter estimate. Typically,
since EEG/MEG signals contain a large amount of noises or artifacts from outside of the brain,robustness is a desirable property in practice. On the other hand, our model includes no explicit
noise term in the generative model of Eq. (1), which is mainly for computational simplicity as
in standard ICA methods.
Graphical representations of the dependency structure in LCMM are given in Fig. 2. Thee) 6=
sources in LCMM are not independent of each other (Fig. 2(a); see also Fig. 1 (c)), i.e., p(s
Qd
e
p(
s
),
in
contrast
tostandard
complex-valued
BSS/ICAmodels.
Unliketherelatedgenj
j=1
eralizations of ICA modeling energy correlations, the sources are even conditionally dependent(Fig. 2(b)), given the higher-order latent variables (here: state k). This is due to our choice of
the Student-t model instead of a Gaussian model. In the limit of ν → ∞ the pdf is reduced
to a complex Gaussian and the sources become conditionally independent (Fig. 2(c)), but theyremain dependent over the whole data set (with state k marginalized out).

2.3

Parameter estimatione by theAccording to LCMM’s generative model, we obtain the pdf of complex sensor signal x
e to xe:transformation of random variables from se ; W, B, η, ν) = | det W|2
p(x

K
X

e ; diag(bk ), ν),ηk Te (Wx

(6)

k=1

where we explicitly indicate the model parameters in the left-hand side (after the semicolon),and B = (b1 , . . . , bK ) and η = (η1 , η2 , . . . , ηK )T collect the coactivation patterns and the state
probabilities, respectively. Note that in Eq. (6), the determinant is squared because the same
2

Note that other heavy-tailed distributions in the complex elliptical family [43] commonly have the robustness
property and thus could also be used as the source model in LCMM. However, the examination of this option isbeyond the scope of this paper.

6

transformation is required for both real and imaginary parts. The pdf can also be expressed as
e ; A, B, η, ν) =
p(x

K
X

e ; AT diag(bk )A, ν).
ηk Te (x(7)

k=1

The model is thus a constrained form of the mixture of multivariate Student-t distributions,in which the state-dependent scatter matrices are tied with common parameter A. It can bereadily seen from Eq. (7) that the probability density does not change if A (or W) and bk
1
1
are simultaneously replaced by AD− 2 (or D 2 W) and Dbk , respectively, where D is any nonsingular diagonal matrix. This is the well-known scaling ambiguity inherent to BSS/ICA; we fix
the scale by setting every column of A to have unit Euclidean norm.
The parameters of interest in LCMM can be easily estimated by the maximum likelihoodPe (t); W, B, η, ν) with respect to the model parameters when
method, i.e., maximizing N
t=1 ln p(xthe latent variables are marginalized out. Importantly, this does not require any approximationin contrast to other hierarchical BSS models. For simplicity, we fix degrees of freedom ν to a
constant (we set ν = 2 in Sections 3 and 4 below; this choice leads to an infinite variance andstrong robustness) and learn the other parameters {W, B, η} since ν typically has only a small
effect on the final solution (at least if set relatively small for ensuring robustness). The detailed
form of the objective function and its derivatives are given in Appendix A.
We propose to use a quasi-Newton method3 to efficiently optimize the likelihood using reparameterization [41] given byexp(λk )ηk = PK,k0 =1 exp(λk0 )(8)

so that the ηk s automatically satisfy constraints ηk ≥ 0 and K
k=1 ηk = 1. On the other hand, we
don’t constrain the scaling of W or B during the optimization, but rescale them after obtaining
the final solution so that every column of A has a unit norm. Any standard optimization software
can be used for solving unconstrained optimization on new parameter set {W, B, λ}. We foundthat this quasi-Newton method is more efficient than the well-known expectation-maximization
method (simulations not shown).
After estimating the model parameters, the (real-valued) sources are separated by s(t) =e(t) = Wxe (t) for any x(t) or x
e (t). This is possible because
Wx(t) or by taking the real part of s
the demixing matrix is the same for both the original and the Hilbert-transformed data. The
states are inferred by computing their posterior probability:Pe ; diag(bk ), ν)ηk Te (Wxe) = Pp(k | x.
Kee ; diag(bk ), ν)
k0 =1 ηk0 T (Wx(9)

The maximum a posteriori (MAP) estimate of the state is given by taking the state that maximizes this posterior, which gives the final result of the model-based clustering of the sourcecoactivations.

2.4

Choosing the number of states by BIC

Another important issue in learning mixture models is the choice of the number K of states orclusters. We use the Bayesian information criterion (BIC) to select the best K minimizingb + M ln N,
BIC(K) := −2 ln L
3

(10)

We used a Matlab implementation [42] of the limited-memory BFGS by Mark Schmidt, available at
http://www.di.ens.fr/˜mschmidt/Software/minFunc.html

7

b denotes the maximum of the likelihood obtained numerically as explained above and thewhere L
number of free parameters M in LCMM is given by M = K −1+d2 +Kd−d. The use of BIC for
the model order selection in mixture models has been extensively studied in statistics. There aretheoretical results regarding statistical consistency [28], and BIC often exhibits state-of-the-art
performance, as shown empirically [46].

2.5

Relation to previous two-layer extensions of BSS/ICA

Many previous attempts have been made to extend BSS/ICA based on Eq. (1), particularly todeal with the residual dependency structures between power s2j or magnitudes |sj | of the sources
often observed in ICA results [21, 22, 47, 26, 27, 50, 17]; these works were not necessarily
concerned with EEG/MEG.
Initial developments made fixed prior assumptions about the dependencies of the sources
without estimating any parameters in the source model [21]. However, we are concerned withmodels in which the parameters in the second layer (i.e., connectivity patterns) are estimatedas well [26, 37, 29]. Typically, these models are based on a (generalized) linear model of squared
sources s2 := (s21 , s22 , . . . , s2d )T , as already proposed by [22]:E[s2 (t) | u(t)] = φ(Bu(t)),(11)

where B and u(t) = (u1 (t), u2 (t), . . . , uK (t))T are the second-layer mixing matrix and source
vector, and φ is a strictly monotonic function that is applied element-wise. Given varianceE[s2j | u], each source sj (t) is then usually assumed to follow a certain probability distribution:either Gaussian [22, 47] or non-Gaussian [26]. The real-valued linear mixing (Eq. (1)) finallyproduces observed signals x(t).
The following are the main problems with these previous models for EEG/MEG connectivity
analysis: 1) they do not properly model the envelopes of the oscillatory sources, and 2) the
likelihood is often intractable without resorting to approximations, because the continuous latentvariable u(t) needs to be integrated out. A recent study [7] on the statistical modeling of natural
movies actually addressed the first issue without resolving the second one.
In fact, a close connection between LCMM and Eq. (11) is implied by the well-known factthat Student-t distribution belongs to the Gaussian scale-mixture family [see e.g., 4, for a realcase]. We can equivalently re-formulate our model by assuming that sej (t) is a complex (circular)
Gaussian conditionally on state k and introducing a scaling variable uk (t) ≥ 0 that follows aninverse Gamma distribution. The conditional variance can then be written:
e|2 (t) | uk (t)]k = uk (t)bk ,
E[|s

(12)

e|2 := (|se1 |2 , |se2 |2 , . . . , |sed |2 )T denotes the squared envelopes. Now consider a simplified
where |scomplex-valued counterpart of the previously used Eq. (11) given by
e|2 (t) | u(t)] = Bu(t).
E[|s

(13)

The LCMM in Eq. (12) has essentially the same form, if we can constrain it so that only a single
variable uk takes a non-zero value at a time.
Hence, although closely related, LCMM has notable differences from previous energy-correlatione|2 instead of the (squared) magnitudes s2 ,models. First, it models the (squared) envelopes |sthus properly dealing with oscillatory sources (together with Hilbert transform). Second, themodel is tractable and fast to learn because it has only one discrete latent variable instead ofmultiple continuous ones.

8

2.6

Real-valued variant of LCMM

To separately evaluate the effect of using complex-valued formulation instead of a real-valued
kind, we also examine a real-valued counterpart of LCMM in our simulation study below. The
model can also be seen as a simplification of a previous two-layer BSS/ICA such that only asingle variable uk in Eq. (11) takes a non-zero value at a time, where nonlinearity φ is set to an
identity function, implying that E[s2 (t) | uk (t)]k = uk (t)bk . More specifically, the real-valued
LCMM is given as a constrained form of a mixture of multivariate Student-t distributions andits pdf is given byp(x; A, B, η, κ) =KXηk T (x; AT diag(bk )A, κ),

(14)

k=1

where the Student-t pdf for real vector s ∈ Rd is generally given byd+κ

Γ( d+κ
1 T −1 − 2
2 )T (s; Σ, κ) =
1
+
s Σ s
.
κ
κ
Γ( 2 )(κπ)d/2 |Σ|1/2




(15)

Again, the pdf can be expressed using W = A−1 asp(x; W, B, η, κ) = | det W|

KX

ηk T (Wx; diag(bk ), κ),

(16)

k=1

and the same quasi-Newton optimization is used for estimating {W, B, η} (specifically with
κ = 1 in Sections 3 and 4 below). Finally, the state posterior is also given by
ηk T (Wx; diag(bk ), κ)p(k | x) = PK.
k0 =1 ηk0 T (Wx; diag(bk ), κ)

3

(17)

Simulations on artificial data

We next quantitatively compare the proposed method with existing approaches for the analysisof EEG/MEG data. We start with simulated EEG/MEG data so that the ground truth is
known and can be systematically controlled; we provide a real EEG analysis in Section 4. The
goal of the simulation study below is to validate the two key ideas of the proposed method: thecomplex-valued formulation and the unified estimation principle of the two stages of analysis.

3.1

Methods

The simulated EEG/MEG signals were created as follows. First, we applied a band-pass filter(9.5-10.5 Hz) to ten Gaussian temporally white signals sampled virtually at 75 Hz to simulate thealpha-range activities. Then these oscillatory signals were jointly amplitude-modulated blockwise in every 2-second window (150 samples) to show the state-dependent coactivation patterns.
For this purpose, we first created vectors bk whose entries were independently sampled from astandard Gaussian distribution but set to zero if negative; this was repeated until
at least twopentries satisfied bjk ≥ 0.05. Then the j-th oscillatory signal was multiplied by bjk with state
k randomly chosen for each block with uniform probability (k = 1, 2, . . . , 5). The sources with
non-zero bjk ’s were actually coactivated, while very small activity levels were avoided with thisprocedure.
So that the amplitudes of these coactivated sources have non-zero (positive) correlations
within each state, they were further modulated globally by a noisy sinusoidal signal, generated
9

States
2

1

3

1

4

5

1

1

0

2

0

4
5

1

6

0

7

1

8

0

9

1

Clean
SNR=0 dB

100State 2

1

3

SourcesState 1

4

5

10Time [sec]150State 3

State 4

State 5

1 2 3 4 5 6 7 8 9 10
Sources(b) Coactivation patterns(a) Simulated sources

Figure 3: Simulated data: Example of simulated sources and states. (a) Ten simulated sourcessj (t) before (black) and after (gray) adding Gaussian noise with signal-to-noise ratio (SNR) of
0 dB. Clean sources (i.e., before adding noise) exhibit one of five different coactivation patterns
bk in each block (separated by vertical dotted lines) with their amplitudes co-modulated withineach block (see text for more details). Top horizontal bars and numbers indicate five states
k = 1, 2, . . . , 5, corresponding to five coactivation patterns. (b) Five coactivation patterns. k-th
panel shows values of bjk for ten sources j = 1, 2, . . . , 10.

10

by sampling from a Gamma distribution Gamma(2ξ(t), 2) 4 with ξ(t) = 0.9 sin(2πf t/75 + φ) + 1
where f = 1 Hz and phase φ was randomly selected. Then Gaussian white noise was added,
where the noise variance was set to have a given value of a signal-to-noise ratio (SNR), definedas the ratio of the variance. Fig. 3 illustrates an example of ten sources before and after adding
the noise. Finally, sources sj (t) were linearly mixed into the same number (i.e., 10) of sensor
signals xj (t) with square mixing matrix A generated randomly from the standard Gaussian
distribution.
The clustering and source separation performances on these simulated data were compared among the following methods: 1) LCMM, 2) LCMM (real), 3) ICA+MixT, and
4) ICA+Kmeans. The first one, LCMM, is the proposed approach, based on the joint maximum likelihood estimation of the two layers of the complex-valued LCMM. The second one,
LCMM (real), denotes the real-valued counterpart of LCMM (see Section 2.6), which canbe seen as laying between our proposed method and the previous two-layer BSS/ICA models.
The latter two, ICA+MixT and ICA+Kmeans, perform two-stage analysis. Both first usede(t); thenthe complex-valued FastICA [3] (with real-valued W) for separating complex sources s
ICA+MixT directly learned the mixture of the Student-t model (as in LCMM) on the separated sources, while ICA+Kmeans performed standard k-means clustering on log-amplitudes
ln |sej (t)|, where the mean log-amplitude over the channels was subtracted at every t to compensate for the global modulation.
We used the adjusted mutual information (AMI) [48] and the Amari index [2] as specific performance measures for clustering and source separation, respectively. AMI5 corrects normalized
mutual information (NMI) between true cluster k and estimated cluster kb for chance agreements:
b max{H(k), H(k)}b (0 ≤ NMI ≤ 1) and
AMI = (NMI − NMI)/(1 − NMI) where NMI = I(k, k)/
NMI denotes the expectation of NMI under random permutations of the cluster labels where
the numbers of clusters and cluster members are unchanged (I and H denote the sample mutual
information and marginal entropy). AMI is thus expected to be zero under this random permutation and is upper-bounded by one; the upper bound is achieved only when the two clusterings
are perfectly matched. The Amari index, which is a standard performance measure for linear
BSS problems, is defined by
d
d X
X

d
d
X
X
χjj 0
χjj 0
Amari index =
−1 +
−1 ,
maxk |χjk |maxk χkj 0
j=1 j 0 =1
j=1 j 0 =1







(18)

where χjj 0 denotes the (j, j 0 )-th element of matrix A−1 Atrue with estimated and true mixing
matrices A and Atrue , respectively. This index is nonnegative and equals zero if and only if thetrue mixing matrix is recovered up to the permutation and scaling of the columns.

3.2

Results

Figure 4 quantitatively compares the performances of clustering (panels on the left) and source
separation (panels on the right) by the above four methods with different numbers of clusters
estimated by the model, K = 2, 5, 8, while the number of true states is always 5. Each boxplot
displays the result of 50 runs in each of the different sample sizes N . The SNR of the sources
was specifically set to 20 dB in this figure. In the panels on the right, “ICA” corresponds toboth ICA+MixT and ICA+Kmeans.
As is clearly seen in the left three panels, LCMM achieved the highest AMI (medians) inevery condition, and LCMM (real) consistently showed a lower AMI. In contrast, on the righthand-side panels, these two methods showed very similar Amari indices. These results imply
4
5

The pdf is given by p(x) = ba xa−1 exp(−bx)/Γ(a) if x ∼ Gamma(a, b).
We used the matlab code available at https://sites.google.com/site/vinhnguyenx/softwares.

11

0.6

Amari index

Adjusted mutual information

60

0.8

0.4
0.2
0
1500

40

20

0

3000 6000 12000 24000
Sample size

1500

3000 6000 12000 24000
Sample size

1500

3000 6000 12000 24000
Sample size

60

1
0.8

Amari index

Adjusted mutual information(a) K = 2

0.6
0.4

40

20

0.2
1500

0

3000 6000 12000 24000
Sample size

60

LCMM
LCMM(real)
ICA

0.8

Amari index

Adjusted mutual information(b) K = 5

0.6
0.4
0.2

LCMM
LCMM(real)
ICA+MixT
ICA+Kmeans

1500

40

20

0

3000 6000 12000 24000
Sample size

1500

3000 6000 12000 24000
Sample size(c) K = 8

Figure 4: Simulated data (with SNR=20dB): Performances in clustering (left) and source separation (right) with different settings of number of clusters estimated, (a) K = 2, (b) K = 5,
and (c) K = 8, evaluated respectively by adjusted mutual information (AMI) [48] and Amari
index [2]. Actual number of clusters was always five. AMI is scaled between 0 (completely
random) and 1 (perfect clustering); Amari index becomes zero if estimated A recovers the true
mixing matrix up to the permutation and scaling of the columns. Thus, on the left, high values
are good; on the right, low values are good. Each panel displays boxplots at different sample
sizes N . Each boxplot indicates median, interquartile range, and entire data range of 50 runs,
excluding outliers indicated by ‘x’. See text for legends of methods.
12that the complex-valued formulation of our LCMM is particularly beneficial for obtaining better
clustering without degenerating the source separation. The two-stage methods, ICA+MixT
and ICA+Kmeans, also exhibited lower (i.e., worse) AMI values than those of LCMM (butnot necessarily of LCMM (real)), while they exhibited higher (i.e., worse) Amari indices. The
two-stage methods are less accurate in both clustering and source separation than the proposed(complex-valued) LCMM method.
Although we obtained the best performance with K = 5 (true number of clusters), therelative performance among the four methods was qualitatively similar for different Ks. Thatis, the proposed method outperforms other methods even when the number of clusters K is
misspecified.
To further examine how the result changes with different noise levels, we also conductedsimulations with different SNRs for generating the source signals. The number of clusters K in
the model was simply set at the true one (K = 5). Fig. 5 shows the result in the same format asthat of Fig. 4. The relative performance of the four methods was qualitatively similar to Fig. 4in every SNR setting. This showed that LCMM improved clustering without degenerating thesource separation over the other methods even when the SNR is relatively low.
Finally, we demonstrated the use of BIC for selecting the number of clusters. Here, we
computed the BIC for K = 2, 3, . . . , 10 and chose the K that minimizes the value. Note that in
the simulation setting here, LCMM does not completely match the true data-generating modeldue to the additional Gaussian noise in the sources. The number of clusters selected thus often
exceeded five, as shown in Fig. 6, while higher SNRs (e.g., 20 or 30 dB) resulted in values closer
to five. In practice, since EEG/MEG usually has a low SNR, these results indicate that the
number of clusters will likely be overestimated by BIC. However, spatial topographies aj also
learned by LCMM can be used to identify and discard such irrelevant clusters that only containnoise or artifacts instead of physiologically meaningful patterns.

4

Experiments on resting-state EEG data

Next, we demonstrate the advantages of using LCMM compared to two-stage methods in areal EEG data analysis. The target data are resting-state EEGs acquired before and after a
BCI-related task, which we expect to contain task-relevant brain states possibly due to mental
rehearsal or retrieval. We examined the patterns (states) found in the resting-state EEGs basedon the labeled EEG data during task as well as their spatial topographies on sensor channels.

4.1

EEG data

Five healthy subjects (three males, two females, 28+/-11 years old) participated in our EEG
experiment. We placed a headcap with EEG electrodes on their heads with electric-conductive
gel SIGNAGEL (Parker Laboratories Inc., Fairfield, NJ, USA) to reduce the impedance of theelectrodes. We positioned 64-channel active electrodes based on the international 10-20 system,
and connected them to an ActiveTwo amplifier (BioSemi, Amsterdam, The Netherlands). An
experimental protocol of this study was approved by the ethical committee at ATR.
The brain activity was measured in each subject during resting states with eyes open andwhile the subjects performed a cued motor imagery/non-imagery task. The experiment consistedof two resting-state (RS) sessions and six task sessions between the two RS sessions. In each
RS session (5 minutes), the subject was instructed to relax without thinking of anything inparticular and without sleeping and to focus on a fixation point at the screen’s center. In eachtask session, the subjects performed a number of task trials in each of which they randomly took
one of the following three actions for three seconds after a visually-cued onset: 1) left: covert

13

0.5

Amari index

Adjusted mutual information

60

0.6

0.4
0.3
0.2

40

20

0.1
1500

0

3000 6000 12000 24000
Sample size

1500

3000 6000 12000 24000
Sample size

1500

3000 6000 12000 24000
Sample size

50

0.8

Amari index

Adjusted mutual information(a) SNR = 0 dB

0.6
0.4

40
30
20
10

0.2
1500

0

3000 6000 12000 24000
Sample size

60

1
0.8

Amari index

Adjusted mutual information(b) SNR = 10 dB

0.6
0.4

LCMM
LCMM(real)
ICA+MixT
ICA+Kmeans

0.2
1500

40

20

0

3000 6000 12000 24000
Sample size

LCMM
LCMM(real)
ICA

1500

3000 6000 12000 24000
Sample size(c) SNR = 30 dB

Figure 5: Simulated data (with K = 5 (true)): clustering (left) and source separation (right)
performance with different noise levels. See caption of Fig. 4 for details.

14

N=3000

N=6000N=12000

N=24000

2 3 4 5 6 7 8 9 10

2 3 4 5 6 7 8 9 10

2 3 4 5 6 7 8 9 10

2 3 4 5 6 7 8 9 10

2 3 4 5 6 7 8 9 10

30 dB

20 dB

10 dB

0 dB

N=1500Figure 6: Simulated data: normalized counts of number of clusters K selected by minimizing
BIC over K = 1, 2, . . . , 10. Rows and columns correspond respectively to different noise levels(signal-to-noise ratios) of simulated sources and to different sample sizes N . True number ofclusters was five.
imagination of a left-hand movement, 2) right: covert imagination of a right-hand movement,
and 3) idle: no imagination of hand movements.
The EEG data were acquired at a sampling frequency of 256 Hz, band-pass filtered off-line
to 1-50 Hz (fourth-order Butterworth, zero-phase), and re-referenced to the common average.
Part of the RS data and some trials in the task data were rejected due to gross contamination.
Typical ocular, cardiac, and muscular artifacts were also identified and removed by FastICA [20]with visual inspection and frequency analysis, which was done separately for each subject andalso for the RS and task data. The data were further band-pass filtered in certain frequency
bands of interest and then Hilbert-transformed. We focused on two frequency bands of interest,8-12 Hz (alpha) and 13-30 Hz (beta), in line with previous neurophysiological studies on motor
imagery [e.g., 38] as well as those on resting-state brain networks [e.g., 6].

4.2

MethodAll the LCMM parameters were estimated from the RS data alone where the two (pre- and post) RS sessions were combined. We emphasize that the task data and labels were not used forthe parameter estimation but only to validate the learned model in a post-hoc manner. Before
the parameter estimation, the RS data were spatially prewhitened and dimensionality-reducedby PCA, so that 99% of the sample variance was kept. Note that the effective dimensionality ofthe data had already been reduced above by removing the artifactual dimensions using ICA. As
a result, the number of sources d was selected as 12, 12, 11, 7, and 31 for the five subjects in thealpha band and 17, 12, 17, 14, and 31 in the beta band. For comparison, we also applied the twostage method ICA+MixT, as explained in Section 3.1, to the same data to examine the effectof unifying the two stages of the analysis by LCMM. In both methods, we ran the algorithm
ten times, each from different initial parameters to converge, for every preselected number ofstates K = 10, 20, . . . , 100. The best K in LCMM was chosen so that the median of the ten
BIC values achieved the minimum, and the same number K was also used in ICA+MixT.
We then evaluated how well each coactivation pattern bk found in the RS data discriminates
the two physiologically different brain states corresponding to the motor imagery (i.e., left andright) and non-imagery (i.e., idle) labeled in the task EEG data. We used a standard performance
criterion for binary discrimination, the area under the receiver operating characteristic (ROC)
15

Fpz Fp2
Fp1 AFz
AF4AF8
AF7AF3F4 F6 F8
F2
Fz
F3 F1
F5
F7
FCz FC2 FC4FC6
FC3 FC1
FT8
FC5
C2
Cz
FT7
C4
C6
C3 C1
T8
C5
T7
CPz CP2CP4CP6
CP1
CP3
CP5
P2 P4 P6 TP8TP7P3 P1 Pz
P8
P5P7
PO3 POz PO4PO8 P10
PO7
O1 Oz O2
P9
Iz

Figure 7: Real EEG analysis: Channel names and spatial layout of 64 electrodes. Layout
corresponds to spatial topographies on scalp shown in Figs. 9 and 10.
e | k) to detect
curve or AUC and evaluated them as follows. We used the k-th signal model p(xwhether a task trial was in the k-th state. This state was supposed to be detected if log-likelihoode | k), averaged over the imagery/non-imagery period after the cued onset (where the initialln p(x
0.5 seconds of this 3-second period were discarded to avoid the transient effect), was above or
below a certain threshold. The log-likelihood of the k-th state, up to the irrelevant scaling andP
ej |2 ), as derived byadditive constants, can be evaluated at each time point by ln(1 + dj=1 b−1
jk |ssetting ν = 2 in Eq. (5) and removing the time-invariant constant terms from its logarithm.
The occurrence of the k-th state may be associated arbitrarily with motor imagery or with
non-imagery for each case of which the ROC curve was drawn by plotting the true positive rate
against the false positive rate at many different threshold values. Thus we had two AUC values(computed by a trapezoidal rule) for each state from which the greater one was simply chosen;
the AUC becomes close to one if the state occurrence discriminates well between motor imagery
and non-imagery, and it is close to 0.5 if the occurrence does not discriminate it at all. Notethat this AUC does not distinguish between left and right because that seemed too difficult inthese data based on our preliminary analysis (not shown here).
To compare these states that exhibited high AUC values, we further examined the timecourses of their log-likelihood and the sensor-level topographies corresponding to them. The
topographies were drawn by evaluating how the occurrence of each state changed the powerin the frequency band of interest at each site of the electrodes. We first obtained a robust
2 ∝ P b a2
estimate of the state-conditional variance (power) for each sensor channel by σck
j jk cj(c = 1, 2, . . . , 64), where acj denotes the estimated mixing coefficient from the j-th source to
P 2
the c-th sensor channel and total variance c σck
was simply normalized to one because it
was undefined due to the choice of ν = 2. The topographies were then plotted by spatially2 − σ 2 )/σ 2 from the “grand variance” σ 2 overinterpolating the percentage deviations 100 × (σckccc
64 electrodes (see Fig. 7 for the layout and channel names). Grand variance σc2 was set to the
2 during the task period, i.e., σ 2 = P η σ 2 , whereexpectation of state-conditional variances σckc
k k ckηk was re-estimated by the posterior state probabilities during the imagery/non-imagery periodsof the task trials.

4.3

Results

The discriminability of each state between the motor imagery and idling trials was examined interms of the sample distribution of their AUC values. Fig. 8 shows the normalized histogramsof the AUC collected from the ten different runs of the parameter estimation, where each panel
shows the two normalized histograms by LCMM (“unified”) and ICA+MixT (“2-stage”) for16

a particular subject and frequency band. The number of samples (AUC values) summarized ineach panel was thus 10 × K, where number of states K was selected by BIC as K = 70, 70, 90, 20
and 90 for the five subjects in the alpha band, and K = 30, 10, 20, 10 and 70 in the beta band.
The figure shows that in most cases, the unified method exhibited greater variation in AUC
than the two-stage method. For example, the range of AUC clearly expanded to both sides ofthe distribution in Subjects 1 (beta), 3 (beta), and 5 (alpha), while mostly to the right (greater)side in Subjects 3 (alpha), 4 (alpha, beta) and 5 (beta). The greater variation in the AUC
distribution, even if it was both-sided, implies that it has a greater chance to contain brain
states that exhibit higher AUC values. The two vertical lines in each panel further indicate
the upper quartiles of AUC for the unified (solid line) and two-stage methods (dashed line),showing consistent increases of the upper quartiles by LCMM from those by ICA+MixT. In
other words, the lower bound of the top 25% of the AUC values consistently shifted to the rightas a consequence of unifying the two stages of analysis.
To get insights into the difference of these high-AUC states between the unified and twostage methods, we further analyzed the average dynamics of the states that exhibited the best(highest) AUC in each subject and frequency band. Figs. 9 and 10 display the trial-averaged
time courses of the log-likelihood of those states obtained for the alpha and the beta bands,
respectively, around the 3-second period of motor imagery/non-imagery. For the alpha band
(Fig. 9), the dynamics clearly differ between the motor imagery (left and right) and the nonimagery (idle) conditions in LCMM (left column), especially in Subjects 1, 2, and 5, withsmaller differences among them in ICA+MixT (right column). For the beta band (Fig. 10), thedynamics again exhibited similar differences between the imagery and non-imagery conditions
in every subject. The results by LCMM (left column) and ICA+MixT (right column) arevery similar in Subjects 2, 3, and 4, but in Subjects 1 and 5 the dynamics again clearly differ
between the imagery and non-imagery conditions in LCMM with fewer differences among themin ICA+MixT.
The corresponding spatial topographies of the power deviations, given in Figs. 9 and 10,
provide further neurophysiological insights in conjunction with those of the averaged dynamics.
For example, Fig. 9 (left column) suggests that the best-AUC states of Subjects 1 and 2 are both
associated with the event-related decrease of frontal alpha power in the idle condition but notin the two motor imagery conditions. This is readily seen in the decrease of the log-likelihood
during the idle condition (green line) from the baseline level and with those topographies thatexhibit positive values on the central areas. The figure also suggests that the best-AUC state of
Subject 5 is associated with the event-related decrease of the alpha powers during the left andright conditions at the bilateral Rolandic (central) areas, which are the major regions of interest
related to motor imagery [38]. This is seen in the log-likelihood increase during motor imagery
in conjunction with the negative deviations in the alpha power around those areas, as indicatedby the topographies. The topographies obtained for the beta band (Fig. 9) seem more difficult
to interpret. Further neurophysiological interpretation is beyond the scope of this paper.
These topographic changes in power, seen at the sensor level, are actually caused by different coactivation patterns at the level of the underlying sources. Fig. 11 shows an example ofthe coactivation patterns (bk ) and the topographies (aj and their element-wise squares) corresponding to each source obtained by LCMM. The top row shows the coactivation pattern of thetwelve sources, obtained for Subject 2 in the alpha band, which corresponds to the best-AUC
state shown in Fig. 9 (left column, Subject 2). Here in this state, sources 3 and 10 have thelargest powers, followed by 5, 6, and 7, and the rest have relatively small powers. With thecoactivation patterns for other states, perhaps these high-AUC states can be characterized by
the relative deactivation of sources 1 and 2 or the relative activation of sources 5, 6, and 7, forexample. Even though we omit the details, the topographies given at the bottom will be useful17Figure 8: Real EEG analysis: Distributional difference of AUC values evaluated for each obtained
state between unified and two-stage methods. Each panel shows two histograms corresponding
to LCMM (unified) and ICA+MixT (2-stage) for particular subjects and frequency bands, as
indicated on left and right sides of figure, respectively. In each panel, two normalized histograms
are superimposed with their overlap shown by transparency. Vertical axis is scaled in each panel,and horizontal axis denotes AUC. Two vertical lines indicate upper quartiles for two methods:
LCMM (solid line) and ICA+MixT (dashed line).

18

Figure 9: Real EEG analysis: Temporal dynamics of log-likelihood of states that exhibitedbest (highest) AUC in Fig. 8 for alpha band (8-12 Hz). Each panel displays time courses of loglikelihood, trial-averaged in each of three task conditions, left (blue), right (red), and idle (green),
for a particular subject and by either LCMM (unified) or ICA+MixT (2-stage), as indicatedon left and top of figure. Solid lines indicate moving-averages using time windows of 0.5 seconds,
where shaded intervals indicate standard deviation in each moving window. Corresponding scalp
topographies shows percentage deviations in frequency-band power at each electrode (Fig. 7 forchannel names) interpolated spatially, which represents how the occurrence of state changes thepower from the grand mean at the sensor level.

19

Figure 10: Real EEG analysis: Temporal dynamics of log-likelihood of states that exhibited best(highest) AUC in Fig. 8 for beta band (13-30 Hz). See Fig. 9 for more details. Moving-average
was taken by time windows of 0.25 seconds.

20

Figure 11: Real EEG analysis: Examples of coactivation patterns bk of sources and spatial
topographies corresponding to each source obtained by LCMM (Subject 2, alpha band). The
ten barplots display coactivations associated with ten different states k. Five from top arethose of highest AUC values, and the rest are lowest AUC values, as indicated at right of eachrow. Vertical length of the j-th bar (from left) in the k-th row represents relative value of bjk
with vertical axis scaled separately in each row. Scalp topographies at bottom show values of
mixing coefficients acj (lower) and its squared a2cj (upper), where blue and red colors correspondto negative and positive signs, respectively. Note that the signs of coefficients acj may bearbitrarily flipped for each j due to indeterminacy inherent to LCMM (or ICA). Numbers j atbottom indicate twelve sources.

21

for further interpretation.

5

Discussion

We presented a novel unsupervised method for non-stationary functional connectivity analysisof EEG/MEG sources. Our simulation studies confirmed that the proposed unified method
often outperformed the conventional two-stage method, in terms of both source separation andclustering performances. Real EEG data analysis also showed that the unified method finds
coactivity patterns that discriminate well between motor imagery and idling states with a higher
probability than the two-stage method. In addition, the proposed method performs clustering inthe source space to give further neurophysiological insights beyond sensor-space clustering [e.g.,5, 44], as demonstrated in Section 4.3.
Our real EEG data analysis suggested that the log-likelihood values of the states, obtainedby our unsupervised analysis, may provide better higher-order signal features relevant to BCI(e.g., for the onset detection of motor imagery) than those by conventional methods. However,as seen in Section 4.3, LCMM often did not consistently improve the discriminability (AUC) of
every state but rather strengthened their contrasts: the states obtained by LCMM may containboth more discriminative (higher AUC) and more indiscriminative (lower AUC) ones than thoseby the two-stage methods. Hence, we must carefully select the discriminative features (states)
while avoiding the indiscriminative ones to successfully apply our method to BCI. Such a feature
selection could be done, e.g., based on neurophysiological interpretations or using additional
task-based experimental calibration. This remains open for future investigation.
Some recent studies have combined BSS with effective (directional) connectivity analysis toanalyze the causality between neural activities, for example, autoregressive (AR) models [13, 16,
see also [12] to solve the inverse problem rather than BSS], generalized autoregressive conditional
heteroscedasticity (GARCH) models [50], or a structural equation model (SEM) [17]. Our ideaof unifying BSS and EEG/MEG connectivity analysis in a hierarchical statistical model is thus
not completely novel, but in the present study, we focused on a different type of statisticalconnectivity: functional connectivity based on envelope correlations. More importantly, we focused on analyzing non-stationary functional connectivity in terms of underlying patterns/states
and their dynamics. This is a conceptually crucial difference from previous studies that wereconcerned with static connectivity, i.e., those unchanging over time.
The LCMM proposed here was shown to be a reasonable statistical model to perform our
specific analysis on resting-state EEG/MEG signals, but obviously it has some limitations fromthe perspective of generative signal modeling. First, each LCMM state can describe only the
positive correlations in the source envelopes, as is readily seen from the interpretation of Eq. (12),where latent factor uk (t) and coefficients bk are both positive. The states themselves are, on theother hand, correlated negatively since they do not occur simultaneously due to the assumptionof the finite mixture model. Hence, our method cannot analyze any brain states characterized bynegative envelope correlations or states that are positively correlated with each other. Second,
the i.i.d. assumption is too simplistic since brain activity is obviously not independent overtime. More general models corresponding to Eq. (13), possibly with nonlinearly φ in Eq. (11) as
well as some temporal dynamics model on latent variables, will make the model more realistic.
However, the model complexity must be carefully controlled to achieve good predictability andto maintain the model’s tractability.
Another practical limitation of our method is that it can currently handle only a single
frequency band of interest. Although it is very typical to limit an EEG/MEG analysis to acertain frequency band, cross-frequency interactions are sometimes of particular interest. Some
recent studies have actually combined a complex-valued ICA with time-frequency decomposition
22

for spontaneous EEG/MEG analysis with more flexibility on the spectral nature [e.g., 24, 39].
A similar technique can probably be used to extend our method to analyze the data in multiple
frequency bands.
Acknowledgments We gratefully thank Motoaki Kawanabe and Okito Yamashita for theirhelpful comments and discussions and Takayuki Suyama, who provided the opportunity forthis research. This work was supported by a contract with the Ministry of Internal Affairs andCommunications entitled, “Novel and In- novative R&D Making Use of Brain Structures” and by
JSPS KAKENHI 25730155. J.H. was partially supported by Strategic International Research
Cooperative Program, Japan Science and Technology Agency (JST). A.H. was supported by
the Academy of Finland, Centre- of-Excellence in Inverse Problems Research and Centre-ofExcellence Algorithmic Data Analysis.

A

Maximum likelihood estimation of LCMM

In our LCMM, the maximum likelihood estimates of parameters of interest {W, B, η} are givenPby the minimizer of the negative (average) log-likelihood L := (1/N ) N
t=1 L(t), where t-th term,e(t); W, B, η) + const., is given byL(t) := − ln p(xL(t) = − ln

K
X





ηk exp −G dXe (t)|2
|wT
jx
j=1

k=1

bjk


−

d
X



ln bjk  − 2 ln | det W|,

(19)

j=1

where G(u) = (d + ν2 ) ln (1 + 2u/ν), according to Eqs. (5) and (6). We minimize L by a quasiNewton method with respect to W, B, and λ using the reparametrization of Eq. (8). Thisrequires the first derivatives of L(t), given by




X

d
 1
|sej (t)|2 |sei (t)|2 ∂L(t)= qk (t)
−g
,
 bik
∂bik
bjk
b2ik 
j=1


(20)



 d

K
∂L(t)  X
qk (t) X|sej (t)|2 −Te∗j (t) + se∗i (t)xej (t)) − 2Wij,=g(se (t)x

 i
∂wij
b
b
ik
jk
j=1
k=1

(21)

∂L(t)
= qk (t) − ηk ,
∂λk

(22)

−T
where g(u) := G0 (u) = (ν +2d)/ (ν + 2u), asterisk ·∗ denotes a complex conjugate, Wij
denotesthe (i, j)-element of the transposed inverse matrix of W. qk (t) denotes the posterior probabilityof the k-th state, given by



X
d

|sej (t)|2
qk (t) ∝ ηk exp −G
bjk
j=1



−

d
X



ln bjk  ,

(23)

j=1

where ∝ means that the left-hand side is proportional to the right-hand side up to a constantfactor independent of k. Notice that Eq. (23) is equivalent to Eq. (9). To obtain Eq. (21) above,
we used the relation given byo∂|sei |2∂ ne ]2 + Im[w T
e ]2
=
Re[wT
i x
i x
∂wij
∂wij

(24)

e ]) Re[xe ]) Im[xej ] + 2(w Tej ]= 2(wTi Re[xi Im[x(25)

ej ] + 2 Im[sei ] Im[xej ]= 2 Re[sei ] Re[x(26)e∗j + se∗i x
ej ,= sei x(27)
23

where Re[z] and Im[z] denote the real and imaginary parts of complex number z, respectively.

References[1] E. A. Allen et al. Tracking whole-brain connectivity dynamics in the resting state. Cerebral
Cortex, 24(3):663–676, 2014.
[2] S. Amari, A. Cichocki, and H. H. Yang. A new learning algorithm for blind signal separation.
In Advances in Neural Information Processing Systems 8 (NIPS1995), pages 757–763, 1996.
[3] E. Bingham and A. Hyvärinen. A fast fixed-point algorithm for independent component
analysis of complex valued signals. International Journal of Neural Systems, 10(1):1–8,
2000.
[4] C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
[5] J. Britz, D. Van De Ville, and C. M. Michel. BOLD correlates of EEG topography reveal
rapid resting-state network dynamics. NeuroImage, 52:1162–1170, 2010.
[6] M. J. Brookes et al. Investigating the electrophysiological basis of resting state networks
using magnetoencephalography. Proceedings of the National Academy of Sciences of theUnited States of America, 108(40):16783–16788, 2011.
[7] C. F. Cadieu and B. A. Olshausen. Learning intermediate-level representations of form andmotion from natural movies. Neural Computation, 24(4):827–866, 2012.
[8] S. Dähne et al. Finding brain oscillations with power dependencies in neuroimaging data.
NeuroImage, 96:334âĂŞ–348, 2014.
[9] F. de Pasquale et al. Temporal dynamics of spontaneous MEG activity in brain networks. Proceedings of the National Academy of Sciences of the United States of America,
107(13):6040–6045, 2010.
[10] S. M. Doesburg and L. M. Ward. Synchronization between sources: Emerging methods forunderstanding large-scale functional networks in the human brain. In J. L. Perez-Velazquez
and R. Wennberg, editors, Coordinated Activity in the Brain, pages 25–42. Springer, New
York, 2009.
[11] K. J. Friston. Functional and effective connectivity in neuroimaging: A synthesis. Human
Brain Mapping, 2:56–78, 1994.
[12] M. Fukushima, O. Yamashita, T.R. Knösche, and M. Sato. MEG source reconstruction
based on identification of directed source interactions on whole-brain anatomical networks.
NeuroImage, 105:408–427, 2015.
[13] G. Gómez-Herrero, M. Atienza, K. Egiazarian, and J.L. Cantero. Measuring directional
coupling between EEG sources. NeuroImage, 43:497–508, 2008.
[14] M. Grosse-Wentrup. Understanding brain connectivity patterns during motor imagery forbrain-computer interfacing. In Advances in Neural Information Processing Systems 21(NIPS 2008), pages 561–568, 2009.
[15] M. Grosse-Wentrup. What are the causes of performance variation in brain-computerinterfacing? International Journal of Bioelectromagnetism, 13(3):115–116, 2011.
24

[16] S. Haufe, R. Tomioka, G. Nolte, K.-R. Müller, and M. Kawanabe. Modeling sparse connectivity between underlying brain sources for EEG/MEG. IEEE Transactions on Biomedical
Engineering, 57(8):1954–1963, 2010.
[17] J. Hirayama and A. Hyvärinen. Structural equations and divisive normalization for energydependent component analysis. In Advances in Neural Information Processing Systems 24(NIPS2011), pages 1872–1880, 2012.
[18] J. Hirayama, T. Ogawa, and A. Hyvärinen. Simultaneous blind separation and clustering ofcoactivated EEG/MEG sources for analyzing spontaneous brain activity. In Proceedings ofthe 30th Annual International Conference of the IEEE Engineering in Medicine and Biology
Society (EMBC2014), pages 4932–4935, 2014.
[19] N. Hironaga and A. A. Ioannides. Localization of individual area neuronal activity. NeuroImage, 34(4):1519–1534, 2007.
[20] A. Hyvärinen. Fast and robust fixed-point algorithms for independent component analysis.
IEEE Transactions on Neural Networks, 10(3):626–634, 1999.
[21] A. Hyvärinen and P. O. Hoyer. Emergence of phase and shift invariant features by decomposition of natural images into independent feature subspaces. Neural Computation,
12(7):1705–1720, 2000.
[22] A. Hyvärinen, P. O. Hoyer, and M. Inki. Topographic independent component analysis.
Neural Computation, 13(7):1527–1558, 2001.
[23] A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component Analysis. John Wiley &Sons, 2001.
[24] A. Hyvärinen, P. Ramkumar, L. Parkkonen, and R. Hari. Independent component analysisof short-time Fourier transforms for spontaneous EEG/MEG analysis. NeuroImage, 49:257–
271, 2010.
[25] A. Hyvärinen, K. Zhang, S. Shimizu, and P. O. Hoyer. Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Research,
11:1709–1731, 2010.
[26] Y. Karklin and M. S. Lewicki. A hierarchical Bayesian model for learning nonlinear statistical regularities in nonstationary natural signals. Neural Computation, 17:397–423, 2005.
[27] M. Kawanabe and K.-R. Müller. Estimating functions for blind separation when sources
have variance dependencies. Journal of Machine Learning Research, 6:453–482, 2005.
[28] C. Keribin. Consistent estimation of the order of mixture models. Sankhyā: The Indian
Journal of Statistics, Series A, 62(1):49–66, 2000.
[29] U. Köster and A. Hyvärinen. A two-layer model of natural stimuli estimated with scorematching. Neural Computation, 22:2308–2333, 2010.
[30] B. J. Lance, S. E. Kerick, A. J. Ries, K. S. Oie, and K. McDowell. Brain-computer interface
technologies in the coming decades. Proceedings of the IEEE, 100(Centennial-Issue):1585–
1599, 2012.
[31] N. Leonardi et al. Principal components of functional connectivity: A new approach to
study dynamic brain connectivity during rest. NeuroImage, 83:937–950, 2013.
25

[32] X. Liu, C. Chang, and J. H. Duyn. Decomposition of spontaneous brain activity into
distinct fMRI co-activation patterns. Frontiers in Systems Neuroscience, 7:101, 2013.
[33] M. Mahot et al. Asymptotic properties of robust complex covariance matrix estimates.
IEEE Transactions on Signal Processing, 61(13):3348–3356, 2013.
[34] S. Makeig et al. Evolving signal processing for brainâĂŞ-computer interfaces. Proceedingsof the IEEE, 100(Centennial-Issue):1567–1584, 2012.
[35] E. Ollila and V. Koivunen. Robust antenna array processing using M-estimators of pseudocovariance. In Proceedings of the 14th IEEE International Symposium on Personal, Indoor
and Mobile Radio Communications (PIMRC2003), pages 7–10, 2003.
[36] J. Onton and S. Makeig. High-frequency broadband modulations of electroencephalographic
spectra. Frontiers in Neuroscience, 159:99–120, 2009.
[37] S. Osindero, M. Welling, and G. E. Hinton. Topographic product models applied to naturalscene statistics. Neural Computation, 18:381–414, 2006.
[38] G. Pfurtscheller and C. Neuper. Motor imagery activates primary sensorimotor area inhumans. Neuroscience Letters, 239:65–68, 1997.
[39] P. Ramkumar, L. Parkkonen, R. Hari, and A. Hyvärinen. Characterization of neuromagnetic
brain rhythms over time scales of minutes using spatial independent component analysis.
Human Brain Mapping., 33(7):1648–1662, 2012.
[40] P. Ramkumar, L. Parkkonen, and A. Hyvärinen. Group-level spatial independent component analysis of Fourier envelopes of resting-state MEG data. NeuroImage, 86:480–491,
2014.
[41] R. Salakhutdinov, S. T. Roweis, and Z. Ghahramani. Optimization with EM and
expectation-conjugate-gradient. In Proceedings of the 20th International Conference on
Machine Learning (ICML2003), pages 672–679, 2003.
[42] M. Schmidt. minFunc: unconstrained differentiable multivariate optimization in Matlab.
http://www.cs.ubc.ca/˜schmidtm/Software/minFunc.html, 2005.
[43] P. J. Schreier and L. L. Scharf. Statistical Signal Processing of Complex-Valued Data.
Cambridge University Press, 2010.
[44] L. C. Shi and B. L. Lu. Dynamic clustering for vigilance analysis based on EEG. In Proceedings of the 30th Annual International Conference of the IEEE Engineering in Medicine
and Biology Society (EMBC2008), pages 54–57, 2008.
[45] S. M. Smith et al. Temporally-independent functional modes of spontaneous brain activity. Proceedings of the National Academy of Sciences of the United States of America,
109(8):3131–3136, 2012.
[46] R. J. Steele and A. E. Raftery. Performance of Bayesian model selection criteria for Gaussian
mixture models. In Frontiers of Statistical Decision Making and Bayesian Analysis, chapter
4.1, pages 113–130. Springer, 2010.
[47] H. Valpola, M. Harva, and J. Karhunen. Hierarchical models of variance sources. Signal
Processing, 84(2):267–282, 2004.

26

[48] N. X. Vinh, J. Epps, and J. Bailey. Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance. Journal of Machine
Learning Research, 11:2837–2854, 2010.
[49] T. O. Zander and C. Kothe. Towards passive brainâĂŞcomputer interfaces: applying
brainâĂŞcomputer interface technology to humanâĂŞmachine systems in general. Journal of Neural Engineering, 8(2), 2011.
[50] K. Zhang and A. Hyvärinen. Source separation and higher-order causal analysis of MEG
and EEG. In Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence(UAI 2010), pages 709–716, 2010.

27