A Non-Parametric Approach for the Activation Detection of Block Design fMRI Simulated Data Using Self-Organizing Maps and Support Vector Machine.

Functional magnetic resonance imaging (fMRI) is a popular method to probe the functional organization of the brain using hemodynamic responses. In this method, volume images of the entire brain are obtained with a very good spatial resolution and low temporal resolution. However, they always suffer from high dimensionality in the face of classification algorithms. In this work, we combine a support vector machine (SVM) with a self-organizing map (SOM) for having a feature-based classification by using SVM. Then, a linear kernel SVM is used for detecting the active areas. Here, we use SOM for feature extracting and labeling the datasets. SOM has two major advances: (i) it reduces dimension of data sets for having less computational complexity and (ii) it is useful for identifying brain regions with small onset differences in hemodynamic responses. Our non-parametric model is compared with parametric and non-parametric methods. We use simulated fMRI data sets and block design inputs in this paper and consider the contrast to noise ratio (CNR) value equal to 0.6 for simulated datasets. fMRI simulated dataset has contrast 1-4% in active areas. The accuracy of our proposed method is 93.63% and the error rate is 6.37%.

Introduction

Functional magnetic resonance imaging (fMRI) has been used to investigate the function of different regions of the human brain. For analyzing fMRI datasets, many statistical analysis and processing methods have been used. These methods contain model-based analysis of variance, correlation method, principal component analysis (PCA), and independent component analysis (ICA).[1] All these methods have the ability to detect active areas.[1] In addition to these mentioned methods, one of the simplest ways for the activation detection is subtraction method.[2] Another common way for the activation detection is cross-correlation (CC) method that uses similarity of real time series of fMRI with reference to time series for the activation detection, and variety of reference signals have different results.[23] K-means, Fuzzy K-means algorithm,[4] Multiple Signal Classification (MUSIC) method,[5] Markov random field,[67] Maximum likelihood ratio (MLE),[8] and wavelet transform[91011] are applied for the activation detection of fMRI datasets. T-test and f-test,[12] canonical correlation analysis,[13] CC analysis by considering trends,[14] Generalized linear model (GLM),[15] and generalized canonical correlation analysis[16] are examples of parametric methods that use a model for active signals. Clustering algorithms, machine learning, PCA, and ICA are examples of non-parametric methods that have not any prior knowledge about signal model.[17] Among all of mentioned algorithms, non-parametric methods that use clustering and classification algorithms have better performance than other parametric related algorithms, because non-parametric methods do not consider any model for hemodynamic response and input signals. There are many supervised pattern recognition methods for fMRI data classification.[18] These methods contain machine learning, K nearest neighborhood, Gaussian naïve Bayes, and linear discriminative analysis.[1920] Classifiers that are based on support vector machines (SVMs) have better performance than others.[1920] For the first time, this method was used on fMRI datasets by Cox and Savoy[19] There are some important reasons for choosing SVM classifier in this work. First, SVM can create a maximum margin in a high dimensional feature space by using kernel function. Second, SVM has enabled many theoretical bounds on the generalization to estimate its capacity, and finally SVM can select only support vectors help us find the appropriate differences between two groups. fMRI images are high dimensional datasets reducing the accuracy of data analysis in the classification methods. Two common ways that have been ever used in dimension reduction and feature extraction/selection-based data classification are PCA[2122] and ICA.[1] ICA is a modification-based data driven approach for extracting independent sources from fMRI datasets. But linearity and lack of having fine performance in representing the results in fMRI datasets are the major problems of this method.[2324] The main obstacle of using PCA is that it merely lessens second order affiliation between each component, and it is supposed that each component is mutually orthogonal, but we do not have any knowledge on whether the physiological meanings of these components are mutually remains orthogonal or not.[24] ICA has better performance than PCA in reducing dependency, but it is unable to select precise components.[23] If the fMRI time series have high spatial/temporal correlations, ICA will not be able to extract the signals.[23] We use SOM for dimension reduction as a feature extraction/selection before accepting SVM for the activation detection in our work. SOM is a data-driven approach, and artificial neural network is used as a non-parametric model for mapping input data from a high dimension into a set of nodes in a low dimensional map. By choosing a large number of nodes, we are able to increase the certainty of proper fitness and flexibility in mapping of the data to low dimension spaces. This allows us to access informative patterns of neuronal communications between nodes in a low dimension state (often a 2-D map). After reducing dimension and labeling datasets by using SOM, a linear kernel SVM classifier is trained on the selected features. This low dimensional data can be treated as features extracted from the original data. Then SVM applied for the activation detection of test datasets. Our non-parametric method is also compared with parametric and non-parametric methods and verifies better performance. Parametric methods that are used in this paper are such as CC<, hybrid wavelet + SOM, hybrid cross-correlation + SOM (CCSOM), hybrid cross-correlation + K-means, and hybrid wavelet + K-means and then non-parametric methods such as hybrid K-means + SVM, K-means, and SOM. In parametric methods input signals have the main role for regional brain activation detection, the final step is that by choosing P-value or randomization methods and considering the error rate, activation detection has been performed. However, in this paper, we replace clustering methods in the final step of activation detection and compare the results. By doing this, we do not have any concerns for choosing the thresholding value, but accuracy rates are considerable in these cases. In CC method thresholding has low accuracy rate than hybrid correlation methods had (CC + K-means and CC + SOM) but in wavelet method, randomization[17] has higher accuracy than wavelet hybrid methods could do. We do our model on block design fMRI simulated datasets and describe SOM and SVM in details.

Materials and Methods

In this section, self-organizing map (SOM) and SVM methodology are reviewed. Then by combing them and using simulated datasets, we evaluate the results.

Self-organizing map

SOM: It is a model of two-layer artificial neural networks that maps high dimensional input datasets to a set of nodes that arrange in lattice.[2526] SOM has two steps: (i) determination of winner node and (ii) updating weighted vectors associated with winner node and some of neighbor nodes. Before training, weighted vectors associated with winner node of map have been initialized. The training expands over several iterations and is based on competitive learning. In each iteration, vector x = [x1,x2,…,x] ∈ R (n is length of time series) is compared with weighted vector nodes: m = [m1,m2, …, m] ∈ R, i = 1, 2, …, N; N is the total number of nodes to determine the winner node, that is, as same as best matching unit (BMU).[26] BMU is the node such that its weighted node has maximum matching with input vector according to the similarity metric. The common linear metric is Euclidean distance that is as follows:

represents the Euclidean norm. x denotes the inputs, m is the ith weighted node, and m is the BMU weight.[26] After the determination of BMU, the weighted vectors associated with BMU and some of their neighbors are updated on map to make them more similar to their input vector:

m(t + 1) = m(t)+h(t)[x(t) − m(t)]      (2)

where t is the iteration number. h(t) represents the neighborhood kernel which controls neighboring nodes that should be updated. After each iteration, magnitude of updating decreases, and neighborhood kernel takes form of Gaussian function:

where r and r are special coordinates of ith node and winner node, respectively. σ is full-width at half maximum Gaussian kernels that determines the number of neighboring nodes which should be updated, and α indicates learning rate. By increasing learning rate, σ and α decrease.[26] It is common to choose the lattice size equal to N; N is the length of time points of each time series.[27] In this paper, we changed the lattice size from [2 × 2] to [14 × 14] and see the results. We see that the best lattice size for the activation detection is [7 × 7] and has the maximum accuracy.

Support vector machine

SVM: One of the supervised learning algorithms is SVM.[2829] SVM classifier works based on the linear classification of datasets. But by generalization of linear SVM, we can have non-linear SVM classifier. The goal of linear SVM classifier is to choose a linear separator to have maximum safety margin in classes. The linearity assumption on the training data is a very strong assumption. But in real world there are not many data sets that can be separated linearly. In this paper, we use linear SVM for detecting active areas. SVM can easily be extended to the non-linear setting by considering kernel functions. In nonlinear SVM, we introduce a kernel function and define it as . Then a deciding function is chosen as:

We obtain a decision surface in feature space whose complexity depends on the number of dimensions of the feature space. The more complex the linear decision surface in feature space (the larger d) can be written as

Functions of the form φ(x)φ(y) are called kernel functions. If we let k(x̄,ȳ) = φ(x̄)φ(ȳ) be a kernel function, we can write f̂(x̄) as follows:

We are free to choose any kernel functions for our method. In this paper, we choose linear kernel function.

The organization of the paper is as follows: fMRI simulated data sets are discussed in the “Simulated Datasets” section. In the section “Proposed Algorithm”, we present details of the proposed approach. The section “Results and Discussion” describes the results of applying the proposed method to simulated fMRI datasets and the final section “Conclusion” concludes the paper.

Simulated Datasets

We use simulated datasets for comparing our proposed method with other methods. If we do not want to use any fMRI toolbox for simulating datasets, we can create these data sets as previous works do by this way.[1718192021222324252627282930] In block designed inputs, two conditions are considered as the “state of during task” and “state of being idle.” These kinds of inputs have been shown in Figure 1. Note that these inputs are only used for simulation datasets, and we consider that we do not have any prior knowledge about inputs as we do not know in real datasets, because our proposed method have been considered as a non-parametric method. For simulated data sets in classification methods, it is common that a series of data (with known characteristics and homogeneity) are created and divided into two parts, and then based on these datasets, classifiers are trained and are tested. For real datasets, because of do not having any prior knowledge about active areas in brain, different active regions we can divide the volume images of fMRI datasets in two groups. The first half is considered as training data sets. Then active areas were labeled by using any feature-extracted methods (PCA, ICA, SOM, etc.) in training datasets. Of course, there are some regions that the algorithm could not distinguish properly. After this, the next half is considered as test datasets.[18] Another common way is Cox and Savoy[19] works on fMRI data sets. However, using a fixed trained and test dataset will bias the results and decrease the reliability of reports. By using k-fold cross-validation, we can overcome this problem. In this work, we use 5-fold cross-validation for validating our classification and clustering methods. Cross-validation is a statistical method of evaluating and comparing learning algorithms by dividing data into two segments: one is used to learn or train a model, and the other is used to validate the model. Here, we create the simulated dataset for SVM test and train datasets, that is, an image (slice) with size of 64×64×180 by 112 active voxels and 180 time points. Each rest (off) and activity (on) durations is 30 s, and these conditions go on until 180 s. Active areas can be obtained by convolving block stimulus inputs with hemodynamic responses that is considered as follows:[17]

Figure 1

Block designed input patterns in simulated dataset

where d = a, a1 = 6, a2 = 12, b1 = b2 = 0.9 s, and c = 0.35.[1] We also add trends by different shapes as introduced in Afshinpour et al.[30] to time series data. Gaussian noise with variance of 9 and mean of 0 were added to time series. Nonactive areas also have Gaussian distribution. There are several definitions for CNR.[31] However, we choose the definition related in to the domain of hemodynamic response and variance of noise as:

A is the difference between signal peak and baseline. Figure 2 shows this definition.

Figure 2

Illustration of the notation in CNR definitions. Hemodynamic response (activation signal) and noise signal. A defines the amplitude of the activation signal and σ indicates the standard deviation of the noise signal

CNR of simulated fMRI images is 0.6 in simulated datasets. Simulated image is illustrated in Figure 3.

Figure 3

Simulated image that is considered as a slice of brain with special active areas and the activation contrast from left to right are 1, 2, 3, and 4%

Proposed Algorithm

For the determination of active areas based on SVM classifier for dimension reduction and feature extraction, PCA and ICA based methods have been proposed.[32] There are many methods for the activation detection and feature extraction in fMRI data classification.[32] However, the high dimensionality of fMRI images causes computational complexity in classification algorithms. In this work, SOM is used for the dimension reduction and feature extraction of datasets. Then SVM is trained based on these features and used detect the active areas of test datasets. Figure 4 illustrates the block diagram of our proposed algorithm.

Figure 4

Block diagram of our proposed algorithm

In this approach, each voxel has 180 time points and size of each slice is 64 × 64. Therefore, simulated datasets have the size of 64 × 64 × 180. At first, we remove the mean of each time series and then normalize them in training and test datasets. As we mentioned above, by using non-linear metric and changing the lattice size from [2 × 2] to [14 × 14] in SOM, projecting and clustering of 3-D datasets by preserving metric relationship of the inputs have been performed. The number of nodes is variable, and each node has particular weight. These nodes have weighted vectors and 2-D feature vectors that are associated with each of them. A 2-D feature vector is generated by mapping the weight vectors on to a 2-D grid of nodes. The purpose of this topology is that the nodes with highest sample density (as an important feature) should be chosen. After this, simulated dataset will be divided into five folds. Four folds were considered as train datasets and the other remaining dataset was considered as test dataset. Then SVM will be trained by this way. By considering 5-fold and changing the test and train datasets in each fold, we evaluate our algorithm. We do not have any prior knowledge about active areas (as in real datasets). Then by using SOM, we reduce the dimension of data to 2-D map and label the SVM training datasets as active and non-active areas. By considering a bi-class problem and labeling the dimensional reduced data {x,y}, i = 1, …, N, y is −1 or 1 (1 indicates active voxels, −1 indicates non-active voxels). We use these labels in test datasets for the activation detection of brain voxels by using SVM. We compare our method with parametric and non-parametric methods. CC,[17] hybrid cross-correlation + K-means, hybrid cross-correlation + SOM, hybrid wavelet + SOM and hybrid wavelet + K-means are parametric methods that consider a model for active signal. Parametric methods use P-value and randomization[17] for the activation detection of brain different regions. Wavelet transform[17] and correlation methods[131415] are examples of parametric methods that use randomization and P-value respectively. However, choosing the best value for thresholding (P-value and randomization) is the major problem of these methods. In this paper, by using hybrid methods there are no thresholding problems in choosing the best number for threshold value. We see that thresholding in CC method has low accuracy than that in hybrid correlation-clustering methods. But in wavelet method, randomization has high accuracy than hybrid wavelet-clustering methods. In CC and wavelet methods after finding parametric similarity, we can separate the active areas from non-active regions. In reality, these methods are using thresholding for the activation detection.[131415] In this paper, by using hybrid CC methods, we show that these methods have high accuracy than the CC method, and using clustering is better than thresholding the parametric statistics in hybrid CC method. It means that in this case, clustering has good performance than thresholding. However, in wavelet method,[17] randomization has better performance than hybrid wavelet methods that was suggested in this work. As we mentioned above, some parametric and non-parametric methods will be compared in this work. K-means, SOM, and SVM are non-parametric methods for the activation detection of active areas in fMRI simulated datasets. In these algorithms, we do not have any prior knowledge about the hemodynamic function. SVM is a sensitive method in high dimensional and noisy datasets and also by combination of clustering methods with this supervised method (SVM) such that we can label the train datasets. In many cases, K-means does not have good performance when inputs have complicated nature.[33] SOM has excellent performance in clustering methods for analyzing fMRI images.[33] Because of comparing single function and complementary function of unsupervised machine learning in fMRI data analysis, we want to have hybrid methods too; hybrid K-means + SVM and hybrid SOM-SVM. We will subscribe these methods in details as follows:

Parametric methods

CC: The reference signals are created with hemodynamic response and Gaussian noise (variance and mean of noise are 9 and 0, respectively). We named reference signal R and after moving its mean and normalization of R we have Re as follows:[34]

By choosing Y as a time series of each voxel after normalization and removing mean, e is defined as:

e is correlation coefficient that shows similarity of reference and real time series of signals. By Z-fisher transform we have:

It is used to generate a normally distributed z with mean of 0 and variance of 1/N-3 under the null hypothesis.[35] By considering null hypothesis and choosing P-value (in this case α = 0.01), we detect active areas in fMRI simulated datasets.

Hybrid cross-correlation + K-means (CCK): in this method, we calculate relation (11), and then apply K-means algorithm for clustering and activation detection of simulated fMRI datasets. The inputs of K-means algorithm are considered e to detect active areas. We use datasets as same as CC method and accept K-means algorithm as:

We define , n is the number of voxels. By considering C clusters for our datasets that C = {c1, c2, …, c} (k < n), and minimizing objective function:[4]

5-Fold cross-validation has been chosen for this algorithm. We can cluster inputs and detect active regions. Number of clusters is five in each fold for activation detection.

Hybrid CCSOM: To compare the common CC method by hybrid CC method by calculating , as mentioned in hybrid CCK we apply SOM algorithm for clustering our datasets. As we can see in relation (11) input datasets are vectors that come from CC method. By using relations (1)–(3) we can detect active areas. In this algorithm the size of lattice changed from [2 × 2] to [14 × 14] in each fold (5-fold), and it can be seen that the best accuracy for activation detection of active areas was achieved at the size of [7 × 7] and the iteration number was considered 200.

Hybrid wavelet + SOM and hybrid wavelet + K-means: We use wavelet transform as used in Hossein-Zadeh et al.[17] Then, instead of using randomization[17] for detecting active areas, we use SOM and then K-means for our simulated datasets. The steps of algorithm are as follows.[17] We consider reference signal R:

R = [R(0)R(1) ... R(N − 1)]      (14)

By knowing that there are low frequency fluctuations that are told trends we choose trend vector:

T = [T(0)T(1) ... T(N − 1)]      (15)

Discrete wavelet transform (db-4) with J levels of resolution (J = ) is applied to reference signal and trends:

R and are jth level wavelet coefficients of the reference and trend signal, respectively.[17] We want to have maximum separation of trends from active signals to have accurate activation detection. Therefore we calculate power of reference signal and trends:

q and b are power of reference signal and trends, respectively. By considering error function:

We choose the j0 levels that minimize E.[17] Then we accept wavelet transform by j0 levels for reference signal and Yi and then remove the means of signals:

After this we calculate statistic that illustrates similarity of the reference and real time series.[17]

Then we use SOM algorithm for the activation detection of simulated datasets. As mentioned in last section, the lattice size was changed from the [2 × 2] to [14 × 14] in each fold and the results by considering iteration time equal to 200 have been shown in Table 1. By comparing this hybrid method with randomization method in wavelet algorithm,[17] it is clear that the randomization has higher accuracy than hybrid wavelet methods. We want to compare SOM performance with K-means algorithm, then we also do K-means algorithm after calculating and by choosing six clusters, detection of active and non-active areas has been done.

Table 1

The accuracy of algorithms by changing lattice size from [2 × 2] to [14 × 14]

Non-parametric methods

SOM and K-means algorithm: clustering algorithms are widely used in fMRI data clustering. K-means and SOM algorithm are the examples of clustering algorithms. As literature shows that SOM has better performance than K-means in fMRI data clustering.[33] By considering data matrix with size of 4096×180, we apply K-means algorithm with number of clusters eight for our data matrix. To compare K-means with SOM algorithm, we accept SOM by changing the lattice size from [2 × 2] to [14 × 14], and the iteration time was considered 200. In these non-parametric methods also 5-fold has been considered.

Hybrid K-means + SVM: In our proposed method, SOM algorithm is used for labeling active areas in SVM training datasets. But in hybrid K-means + SVM, we use K-means algorithm for labeling training datasets and extracting the features. After this, we do SVM algorithm for the activation detection of fMRI datasets. Next section discusses the performance of our proposed non-parametric method by comparing other parametric and non-parametric methods.

Results and Discussion

In this paper, we present a SOM-SVM based method for the activation detection of simulated fMRI time series. Figure 5 shows the results of our proposed method and other methods. As we can see from Figure 5 our proposed method has better performance than other methods.

Figure 5

Activation detection of fMRI simulated datasets with CNR = 0.6 and contrasts 1–4%. (a) Hybrid SOM + SVM. (b) Hybrid K-means + SVM. (c) Hybrid wavelet + SOM. (d) Hybrid wavelet + K-means. (e) Hybrid CC + SOM. (f) Hybrid CC + K-means. (g) SOM. (h) K-means. (i) CC. (j) Randomization wavelet method

All of the SOM related algorithms in Figure 5 have the lattice size of [7 × 7]. Table 2 illustrates the accuracy and Table 3 shows the error rate of methods. Definition of accuracy and error rate can be given as follows:

Table 2

Accuracies and standard deviation of accuracies in algorithms

Table 3

Error rate (%) and standard deviation of errors in proposed algorithm and other algorithms

True positive (TP) and false positive (FP) are the non-active areas that are detected as a non-active and active areas, respectively. True negative (TN) and False negative (FN) are the active areas that are detected as an active and non-active areas, respectively. As we can see from Tables 2 and 3, our proposed method has the maximum accuracy and minimum error rate among other algorithms. By changing the lattice size of SOM related methods, we can see that the [7 × 7] lattice size has the best performance than other sizes. In Tables 1 and 4, we changed the lattice sizes and see the accuracies, and also the standard deviation (SD) of accuracies of methods has been proposed, respectively.

Table 4

The standard deviation of accuracies in algorithms by changing lattice size from [2 × 2] to [14 × 14]

As we can see, when the lattice size has the minimum size, the accuracy of algorithms is very low. By increasing the lattice size, algorithms have high accuracies. However, in previous works, it was common to choose the lattice size equal to length of fMRI time series. In this work, the length of fMRI time series is 180. Then, we expect that lattices by size of 169 ([13 × 13]) or 196 ([14 × 14]) have best performance than others. But as we see from the results in low and very high lattice sizes, the accuracy will decrease, and we see that the best lattice size that can create high accuracy is [7 × 7]. The lattice sizes changed from [2 × 2] to [14 × 14]. As we mentioned above, it is common to choose the lattice size equal to length of each time series of fMRI datasets. Then, as we see from Table 1, for all lattice sizes, the accuracies of our proposed method are higher than other methods. However, in comparison with other methods, in Table 1 for some lattice sizes, for example, the lattice size of [12 × 12] in SOM + CC has high accuracy than SOM + SVM in lattice size of [2 × 2], or wavelet + SOM method in lattice size of [11 × 11] has high accuracy than the SOM + SVM in lattice size of [2 × 2]. However, we must consider that all the methods should be compared in same lattice sizes. For the comparison of the results in real fMRI datasets, we must choose the accurate lattice size (length of each time series), and it is obvious that in high lattice size ([13 × 13]), that is almost as same as length of time series (180), our proposed method has high accuracy than other methods have. In real data sets, because of choosing a high lattice size, it is obvious that the accuracy of SOM + SVM based method will be higher than other methods and for having the best choice in accurate lattice size in real datasets, we should use the growing SOM, and of course, this selected size is close to the length of time series.

Conclusion

Parametric and non-parametric methods were widely used in fMRI data analysis. Pattern recognition methods are very important in clustering and classification of brain activities. Model-based methods that use similarity between any two signals such as CC, GLM, canonical correlation, and any generalized correlation methods are commonly used, and other methods, which do not consider any model for active signals, have important role in brain data analysis. In this paper, we work out both parametric and non-parametric methods and compare them. Our proposed model has an accuracy of 93.63% and better performance than other compared parametric and non-parametric methods, and also we have the lattice size of [7 × 7] and running time of 5.18 s. For real data sets, we can do by this way or use fMRI images of two persons with different tasks, for example, auditory and visual tasks. Then, one of them was considered as a test input and another as a train by using of nonlinear statistical relationship on SOM for real datasets. In future works, we want to use growing SOM with SVM on real datasets and comparing the results.

Financial support and sponsorship

Nil.

Conflicts of interest

There are no conflicts of interest.