Independent Component Analysis of Functional Magnetic Resonance Imaging (fmri) Data: a Simple Approach

Independent Component Analysis (ICA) separates spatial and temporal components of fMRI data which may consist of activation patterns, cardiac and respiratory tasks and other artifacts. In this study sources of (fMRI) data are separated using ICA based on simple fixed point iteration method and steepest ascent method. Both are the simplest methods used in optimization. However in this study complete matrix W (un-mixing matrix) is updated in each iteration instead of vector based updating of W. This makes the source separation process very fast. Simulated fMRI data is processed using the proposed method and the results are compared with other ICA approaches in terms of speed and accuracy.


INTRODUCTION
Functional fMRI is a tool used for measuring the functionality of the brain by measuring the oxygenation in the blood flow (Pekar, 2006).Usually experiments are performed for fMRI data acquisition, which are based on a block design comprising cycles of activity and non-activity.During such experiments signals from thousands of brain voxels are taken every few seconds, thus forming a complicated time series of mixtures comprising of task activated, non-activated, transiently task activated, respiratory and cardiac functions activated voxels (Martin et al., 1998).Also the intensity of the activated voxels is too low, normally in the range of 10-15% of variance at 1.5T scanner.This low SNR data is normally first de-noised using some special techniques (Amir et al., 2012).Thus extraction of activated voxels from this complicated mixture of data is a very challenging statistical task.
Changes in fMRI signal, including Blood Oxygen Level Dependent (BOLD) are detected using different techniques including subtraction, correlation, time frequency analysis, Principal component analysis (PCA) and ICA etc.
In subtraction or in general correlation based techniques (Bandettini et al., 1993), the voxels which are responsible for activation of the experimental task are extracted using a priori knowledge of the task related model.A correlation of the experimental model and each time course is performed and the spatial contents are termed as active regions in case of strong correlation results.Other approaches including Statistical Parametric Mapping (SPM) (Friston, 1995) uses t-test and f-tests which are uni-variate methods.In time frequency analysis (Mitra et al., 1997), fMRI signals are recorded in frequency domain.This type of analysis is based on the fact that different experimental task related and function related voxels time courses exhibit different frequencies and thus can be distinguished.PCA is another technique used for fMRI analysis (Moeller and Strother, 1991).This technique finds the orthogonal Eigen-images having the greatest variance in the data.Normally the numbers of principal components representing the data with a specific accuracy are less than the actual dimension of the data.Thus PCA can be used for dimension reduction as well.ICA is another technique which is also used for source separation of fMRI as suggested by (Martin et al., 1998;McKeown et al., 1998) and others.ICA is a statistical method which can convert a multidimensional vector into components which are statistically independent (McKeown, 1998), thus making it useful for blind source separation problem and is widely used in other fields as well (Zhiguo et al., 2006;Hong-Bin et al., 2008;Yu and Cheng, 2012).
In this study fMRI data sources are separated using ICA with a matrix based updating rule which makes the updating process very fast.The proposed techniques prove to be efficient in terms of processing time and accuracy on simulated fMRI data.
ICA model: ICA model can be described as: in which X is the observed data of dimension n×m, A is a mixing matrix of dimension n×k and S consist of independent sources of dimension k×m.
ICA problem can now be described as finding the mixing matrix A such that S can be calculated as: and where, W is the un-mixing matrix.W can be obtained by optimizing a cost function with some constraint.
Commonly used cost functions are based on info-max, maximum likelihood, mutual information and kurtosis etc.
For example if kurtosis is used as a cost function (Aapo et al., 2001), then the cost function needs to be maximized with the constraint W T W = I.Since the sources are independent and the kurtosis of non Gaussian sources is non-zero, therefore we need to maximize the kurtosis which is given by: Different flavors of ICA have been proposed in literature; some of them are Fast-ICA (Aapo and Oja, 1997), JADE (Cardoso, 1999), Extended Info-Max (Lee et al., 1999) and RADICAL ICA (Miller and Fisher, 2003).Further details of these algorithms can be found in the referred literature.
Pre-processing of the data: For implementation of most of the ICA algorithms the observed data needs to be centered (zero mean) and whitened.Centering of the observed data is done by subtracting mean from the observed data Z: where, Z the raw is observed data and  ′ is the centered data.
Similarly whitening is done by multiplying the observed data with some whitening matrix V so that the correlation and covariance matrix of Y becomes identity matrix ie E[YY T ] = I.This process is done using the Eigen-decomposition (Aapo et al., 2001) and is given as.where, V = ED -1/2 E T and EDE T = E[ ′  ′ P T ] Another preprocessing step which is specifically done in case of fMRI due to its high dimensional data is the dimension reduction step.Here we have used Singular Value Decomposition (SVD) as a dimension reduction technique (Aapo et al., 2001).

Fourth order contrast function ICA:
In this study we are using to maximize the fourth moment of Y as the cost function as suggested by Aapo et al. (2001): For further simplification the expectation operator E is also omitted and it becomes: The steepest ascent approach to maximize this function with constraint of W T W= I is given by: where, W(n) = The old value of the un-mixing matrix W W(n+1) = The new value of the un-mixing matrix W X = The observed data matrix Another approach is using fixed point iteration method (Aapo, 1999).In the same way as above the method can be implemented on matrices.
To proceed with fixed point iteration method let f(W) = [Y 4 ] = 0. Now adding W on both sides leads to W(n+1) = W(n) + [Y 4 ].But since we are dealing with matrices therefore care should be taken about the dimensions of the matrices.It is evident in the above equation that we are updating the matrix W while the last term is basically the multiplication of the matrix W with the data i.e., X.Since we need that specific W which gives us f(W), therefore we must multiply the last term by X T thus getting the required form of the update equation: It should be noted that +ve sign here indicates that we are maximizing the function.
It should also be noted that W and X are matrices and not vectors and all iteration are performed on matrices which makes the process very fast.Otherwise in case of vectors based iteration we have to iterate vector by vector, which is the case in most of the ICA algorithms and thus making them slow.

SIMULATION AND RESULTS
For the validity of the proposed schemes, simulated fMRI data of Machine Learning and Signal Processing lab (MLSP) is used which is freely available on the web (Correa et al., 2005).Sources and corresponding time courses are shown in Fig. 1, the corresponding data matrix X consist of 100 images each having 3600 pixels.This makes a 100×3600 data matrix X. Dimension of the data is reduced to 8×3600 using SVD while Pre-processing (centering and whitening) of this data is done using Eq. ( 4) and ( 5).This simulated data is first processed using the depletion based ICA (Ruck et al., 1998) and the extracted sources and time courses are shown in Fig. 2. It can be seen that the quality of the sources extracted is not bad, but the quality of time courses is poor.Figure 3 and 4 shows the sources and extracted time courses by Radical (Miller and Fisher, 2003) and MS based ICA (Molgedey and Schuste, 1994) algorithms.Both these methods results are moderate and took execution time 38 seconds and one second respectively.
Finally the simulated fMRI data is processed using Fast ICA (Aapo and Oja, 1997).
Proposed steepest ascent (matrix based) and Fixed point iteration (matrix based) algorithms respectively.Results are shown in Fig. 5, 6 and 7 having the execution time of 5, 0.8 and 0.6 sec, respectively.Figure 8a and 8b shows the kurtosis vs. iterations and error tolerance vs iterations respectively of the steepest ascent (matrix based) method, while Fig. 9a and 9b shows the kurtosis vs. iterations and error tolerance vs. iterations respectively of Fixed point (matrix based) method.Execution time and correlation analysis of extracted sources and time courses of all the methods are presented in Table 1, which shows that average performance of both steepest ascent (matrix based) and Fixed point (matrix based) is better than ICA methods which are compared in this study.

CONCLUSION
In this study we have presented simple ICA methods (based on steepest ascent and fixed point iteration matrix Based) for the analysis of fMRI data.
These methods are simple and can be implemented on high dimensional data since a dimension reduction step is used before its implementation.Also its convergence time is very less due to its mathematical simplicity and its matrix based updating rules as compared to other conventional ICA methods.Quality performance of these algorithms is also comparatively good as can be seen from Table 1.These method have been tested on simulated fMRI data, however they can be tested on other similar problems.

Fig. 1 :Fig. 2 :
Fig. 1: Sources and corresponding time courses Fig. 8: (a) Kurtosis of sources vs iteration and (b) Error tolerance vs iteration (steepest ascent matrix based) Fig. 9: (a) Kurtosis of sources vs iteration and (b) error tolerance vs iteration (fixed point matrix based)

Table 1 :
Execution time and correlations of (S/T) Sources/Time courses with the actual (S/T) Sources/Time courses