Dominant Component Analysis of Electrophysiological Connectivity Networks

2.1 Projective Non-negative Component Analysis (PNCA)

We hypothesize that each connectivity matrix obtained from recording the brain activity of a subject is a linear combination of several fundamental connectivity matrices called connectivity components. Due to the symmetry of connectivity matrices, a vector of all elements of the upper triangular part of any connectivity matrix is considered as a representative of that matrix, and is used as an observation vector y_i for the corresponding subject i. To compute the connectivity components whose mixture approximately constructs the observed connectivity matrices, a matrix factorization model is used as follows.

where columns of Y_p×q, i.e. y_i (1 ≤ i ≤ q), are the SL connectivity matrix representatives, and columns of W_p×r, i.e. w_j (1 ≤ j ≤ r), are representative of the basis connectivity components, i.e. the upper triangular elements of the matrix of the corresponding connectivity component. These components (w_j) are then mixed by the elements of each column of the loading matrix Φ_r×q to approximate the corresponding column of Y [10, 8].

Assuming that Φ is the projection of Y onto W, i.e. Φ = W^TY, the nonnegativity constraint on the elements of W and Φ makes our non-negative component analysis an optimization problem of minimizing the cost function F(W) = ||Y − WW^TY||² with respect to W, where ||.|| represents the matrix norm. Considering the Frobenius norm, the optimization problem can be denoted by [10]

The above cost function can be minimized by a gradient descent approach, i.e. updating $W_{i,j}=W_{i,j}-{\eta }_{i,j}\frac{\partial F}{\partial W_{i,j}}$ with a positive step-size η_i,j where

The gradient descent updating, as stated above, does not guarantee to keep the elements of W non-negative. However, due to the positivity of the elements of Y and by applying positive initialization to W, our non-negativity constraint is guaranteed by having the step-size as follows

which results in the following multiplicative updating procedure [10]

For stability of the convergence, at each iteration, W is normalized by $\mathit{W}=\frac{\mathit{W}}{\Vert \mathit{W}\Vert }$. Starting by initial random positive elements on W, the iterative procedure will converge to the desired W ≥ 0 whose columns are PNCA. Each obtained component w_j (the j^th column of W) is then normalized by its norm ||w_j||, and correspondingly the coefficients, i.e. elements of the matrix Φ = [φ_ji], are adjusted as ${\phi }_{ji}=\frac{{\phi }_{ji}}{\Vert {\mathit{w}}_j\Vert }:1<i<q$.

The resulting non-negative coefficients are an indicator of the weight of the corresponding component in reconstructing the original connectivity matrices. Therefore, we rank each component w_j based on the average of corresponding coefficients, i.e. $\frac{1}{q}{\sum }_{i=1}^q{\phi }_{ji}$.

2.2 Group PNCA Model for SL Networks

As stated by (1), the q connectivity observations, i.e. y_i : 1 ≤ i ≤ q, in the matrix Y are approximated by

Each observation vector per subject i is thus, approximately reconstructed by

Let us suppose, with no loss of generality, that the first q₁ elements are from the first group (e.g. population of ASD) and the remaining from the second group (e.g. TD group). Thereby, the role of each component w_j in reconstructing the corresponding connectivity vector in the first group, i.e. y_i : 1 ≤ i ≤ q₁, is characterized by the corresponding coefficients φ_ji; and so forth for the second group. Therefore, the statistical significance between the set of {φ_ji : 1 ≤ i ≤ q₁} and {φ_ji : q₁ + 1 ≤ i ≤ q} describes the importance of the corresponding basis connectivity component w_j in differentiating the two groups.

2.3 Synchronization Likelihood (SL) Connectivity Networks

Time-frequency synchronization likelihood assumes that two signals are highly synchronized if a pattern of one signal repeats itself at certain time instants within a time period while the other signal repeats itself at those same instants [9]. Such signal patterns at a time instant t_i of channel k can be represented by an embedding vector x_k,ti = [x_k(t_i), x_k(t_i+l), …, x_k(t_i+(m−1)l)] where l is the lag and m is the length of the embedding vector. l and m are typically set to $l=\frac{f_s}{3h_f}$ and $m=\frac{3h_f}{l_f}+1$ where f_s is the sampling frequency, and h_f and l_f are the high and low frequency contents of the signal, respectively [9]. At each time instant t_i, the Euclidean distance is measured between the reference embedding vector x_k,ti and the set of all other embedding vectors at times t_j, i.e. x_k,tj, where t_j lies in the range $t_i-\frac{t_{w2}}{2}<t_j<t_i-\frac{t_{w1}}{2}$ or $t_i+\frac{t_{w1}}{2}<t_j<t_i+\frac{t_{w2}}{2}$ (in this work t_w2 was set to 10 sec and $t_{w_1}=\frac{2l(m-1)}{f_s}$). Then, n_ref nearest embedding vectors x_k,tj are retained [9].

This procedure is conducted for each channel k and each time instants t_i. Then, the SL between channel k₁ and channel k₂ at time instant t_i is the number of simultaneous recurrences in the two channels divided by the total number of recurrences, i.e. SL_ti = n_k1k2/n_ref. The synchronization likelihood at all the time instants t_i are then averaged yielding the final SL between the two channels.

3.1 Simulation Experiments

In order to demonstrate the effectiveness of the proposed method of PNCA in identifying the true underlying connectivity components, we apply our method to simulated connectivity matrices. We compare this with ICA which is the most widely-adopted method for similar purposes. The simulation is performed by using random non-negative numbers as the elements of three 10 × 10 simulated symmetric SL matrices plus a small background random (non-negative) noise. Ten linear mixtures of the simulated connectivity matrices are composed by a random mixing matrix with non-negative elements. The upper triangular part (excluding the diagonal) of the 10 connectivity matrices are vectorized to form the 10 columns of Y_45×10. We apply the fast ICA algorithm [11] as well as the proposed technique in Sect. 2.1 to solve for r = 3 normalized components as columns of W_45×3. For visualization, the SL matrices and the elements of the resulting components from fast ICA and PNCA algorithms are displayed by grayscale images. The results are shown in Fig. 1. It can be seen that while our method produces the right components (as compared to the original components in Fig. 1(a)), ICA is unable to resolve the components correctly. This is due to the fact that ICA forces the components to be statistically independent (i.e. with covariance of identity) whereas PNCA forces the components to be non-negative and orthogonal which leads to localized components due to sparsity. Please note that the resulting components are not ranked and ordered in Fig. 1.

3.2 PNCA on MEG SL Connectivity Networks