Visual Exploration of Functional MRI Data

1. Introduction

New brain imaging techniques, such as functional Magnetic Resonance Imagery (fMRI), allow for recording and analysis of brain activity over time. A fMRI image is a great source of information about brain behaviour; it is a considerable amount of data (approximately 300000 voxels, "three-dimensional pixels", for which between 100 and 1000 observations are collected). More information about fMRI concepts, studies, data and applications can be found on a publicly accessible repository [1].

From the view point of the data mining, the brain is the most complex object to analyze. In general, the identification of the voxels of the brain that represent real activity is very difficult because of a weak signal-to-noise ratio and of the presence of artefacts. The first tests of the current classification algorithms in this field showed that their performances and their qualities of recognition are weak [2]. Because of the difficulty caused by a very large amount of registered data, the main stream of the research projects is focused on testing a model of brain behaviour by the means of univaried statistics. This is a principle of image processing software such as Statistical Parametric Mapping (SPM) [3], AFNI [4] or BrainVoyager [5], which consists of highlighting the more active voxels under comparable conditions. The statistical methods are powerful, but cannot provide conclusions apart from those prefixed by the model. Using these methods, results must inevitably be anticipated, which is not always possible.

In this article, a new interactive data-driven approach to fMRI mining will be presented. The concept of data mining appears useful as complement or as replacement of the classical methods when it is difficult to predict what will occur during acquisition. In our system, a number of clustering methods have been implemented within an interactive tool to emphasize the active zones without having to use a model. The originality of the approach is not only due to real-time clustering, but also to the insertion of domain knowledge and interactivity, directly integrating the expert-physician into the process of the discovery of functional zones in the brain and their organization. In general, the fMRI data can be analyzed from various perspectives. The three dimensional data recorded temporally can correspond to one or several patients or to the history of the same patient, and, finally, it can be completed by medical knowledge. This article focuses on the first two aspects, where the spatial and temporal dimensions of brain activity are crucial. Important consideration is given to the high volume of data, the processing time, and the noise of image acquisition.

The article is structured in the following way. In section 2, the data mining system architecture oriented toward brain exploration is described. Section 3 describes the three implemented clustering algorithms: Kohonen’s Self Organizing Map [6], Fritzke’s Growing Neural Gaz [7], and Goua’s Clustering using Representatives (CURE) [8]. Section 4 presents the principles of user interface and tools for brain visualization. Section 5 illustrates cases of studies carried out on synthetic and real data. The final section concludes the first experiments and indicates further research perspectives.

2. Data-guided approach

The proposed interactive exploration of fMRI images can be classified as a data-guided approach assisted directly by an expert knowledge and by gathered experience. The process of data mining is composed of five phases: the acquisition and selection of the data, pre-processing, clustering, the extraction of rules and concepts, and validation [9].

The source fMRI data comes from the scanner in a form of sequences of 3D images. The cerebral activity is registered as variations of voxels intensity over time. Typically, the patient is never completely motionless and, moreover, other factors interfere with the signals of interest. Therefore the specific pre-processing must then be adapted to each identified artefact. In the clustering phase, classes are created that are composed of voxels with similar behaviour in time. In the knowledge extraction phase, the classification rules are generated describing each cluster of voxels. Once validated, the rules are saved in the knowledge base, and if required, can be reused in following diagnosis. Figure 1 presents a simplified functional schema of the system, from the fMRI acquisition to interpretation and the validation by the expert-physician, who can interact with the system on all phases.

In many cases the interpretation of fMRI data has to be done very quickly. Frequently, the physician, after obtaining the preliminary results of image clustering, is forced to change the acquisition parameters and redefine geometrical or temporal parameters (resolution, zoom, etc). The clustering and cluster explanation (shown on the right part of fig.1) help the expert to discover and understand the generated classes, and if necessary to modify the experiment directly. These modifications are based on the assumptions about cerebral activities, knowledge of brain anatomy, or about other sources of medical information.

The work presented in this article is primarily focused on the clustering phase, emphasizing interactive and dynamic aspects of algorithms of unsupervised learning. In general, it consists of regrouping voxels that have similar characteristics and behaviours into a limited number of relatively homogeneous clusters. Many clustering algorithms have already been applied to fMRI data. The most common are statistical methods such as K-means [9, 10, 11, 12, 13, 14, 15], Principal Component Analysis [15, 16,], and Independent Component Analysis [18, 19, 20, 21, 22]. Interesting results have been also obtained using fuzzy classification [17], hierarchical classifications [17; 23, 24], Kohonen’s Self-Organizing Maps [25, 26]. The advantage of these methods consists of a higher level of interpretation, but these algorithms are very costly in terms of computing time and memory space.

Figure 1.

Schema of interactive exploration of fMRI data

In our research, we have been focused on real-time clustering algorithms able to discover as quickly as possible classes of voxels in fMRI data, allowing experts to insert their preferences, medical knowledge and spatio-temporal constraints. The interested reader may find more technical details in our previous reports and publications [27, 28, 29, 30].

3. Clustering algorithms

Clustering algorithms usually depend on a distance. A distance between voxels has to be defined. A 3D distance between voxels is irrelevant to identify voxels having the same activity. Taking the 3D distance into account would make close voxels - close from a 3D perspective - look more similar than far-away voxels having the same activity. The distance between voxels should only be defined according to their activities. It should be noted that a clustering based on the activity of the voxels without any influence of their localisation is clearly different from segmentation techniques also used to identify areas in fMRI images that relies on a comparison of neighbouring voxels.

The fMRI data are very highly noisy. The sources are many: heterogeneity of the magnetic field, thermal noises, thermal noises of examined tissues, head movements, eyes movements, breathing, internal movement related to the blood flow pulsation [27]. This list of noise disturbances is only related to the image acquisition, but not associated with the sensory or cognitive noise. Therefore a number of pre-processing operations are required to obtain the data to be analyzed. The most commonly used pre-processing refers to co-registration (correction of movements), rephasing of the brain cuts, normalisation, smoothing, spatial and temporal filtering, segmentation. More about the pre-processing methods can be found in [1, 2, 3, 4, 10].

Before clustering, the representation of each voxel has to be carefully chosen. The attributes to describe the sequence of voxel values have to be selected in such a way that the relationships and distances between voxel activities will be well established. But in clustering, the cause of difficulties may be weak a priori knowledge. The same set of data can be differently clustered depending on selected attributes and distance measure. Note that the activity of a voxel is a continuous signal. Therefore a sampling method for a signal is extremely important. One approach consists in generating different attributes describing the signal, e.g. its average, minimum and maximum values, and then using those attribute-values in a traditional attribute-value clustering system. In such an approach, the built-in distance measure of the clustering system is calculated on the intermediate attributes, e.g. the Euclidian distance between the respective average, minimum and maximum values of each voxel. Its success depends on how well the built-in distance and the generated attributes fit together. Lots of attributes can be tested. For instance, in [3] applied wavelets to transform the signal, though they made use of hidden Markov models rather than a clustering technique.

We considered an alternative approach, where the distance is directly calculated on the fMRI data. The fMRI data are transformed into a time series of voxel intensity variations relative to its average as follows: I^a_ave = 1/n Σ I^a_i where I^a_ave is an average intensity of voxel a of a series of n images; S_a = {δ₁, δ₂,…, δ_n}, δ_i = I^a_ave - I^a_i where S_a is fMRI signal.

The distances between two fMRI signals S_a and S_b may be computed as Euclidian distance: d_E = √(δ^a _i− δ ^b _i )² or Manhattan distance: d_M = Σ |(δ^a _i− δ ^b _i )|.

In the system different clustering algorithms can be easily developed. Currently, the five algorithms are available, notably K-means, LBG, CURE, and two neural models: Kohonen’s Self Organizing Map, and Fritzke’s Growing Neural Gaz. The algorithms K-means and LGB are well known and described in many publications. Algorithms like ICA and PCM separate the fMRI signals into a set of well defined components, but have to deal with constraints of their independency and orthogonality. Therefore in the paper, the description of algorithms is only given to the three less known algorithms that enable the expert interactively to improve his or her understanding of the human brain.

Self-Organizing Maps is a topology-preserving clustering algorithm that maps high-dimensional fMRI data into low-dimensional space [24, 25]. SOM creates the map that represents the fine cluster structures and cluster relations. The specification of SOM algorithm is given below:

Algorithm Self-Organizing Map; Parameters: t: time units; tmax: duration of computing; drs = |i-k| + |j-m| : Manhattan distance between two classes r=aij and s=akm; σi, σf, εi and εf : initial and final coefficients of adaptation; σi = σi .(σi/σI) ^t/tmax : neighbouring coefficient; εi = εi .(εi/εI) ^t/tmax : adaptation coefficient; hrs = exp (-d^2rs / 2σ(t)^2) : neighbouring function between classes r and s; l * h : size of the grid ; number of classes; Procedure; 1.Choose parameters values: size of the grid l * h, the duration tmax, the adaptation coefficients: σi, σf, εi and εf; 2.Grid initialization taking values respecting neighbouring class proximity; 3.Select at random an input signal ξ.; 4.Search for a=g(ξ) of the winning class of ξ; the closest vector of reference; 5.Adapt each class according to the formula wa =wa + ε(t) . haξ . (ξ - wa); 6.Increment the time t=t+1; 7.If t<tmax, then return to step 3, else stop.

However, its fixed topological structure would not help in our application, since there is no a priori topological relationship between the classes. The problem to solve concerned the validity of discovered clusters and choosing the number of selected clusters.

Thus we have preferred the Growing Neural Gas (GNG) algorithm [7]. Its main advantage is that the number of classes is not fixed in advance, as in most clustering algorithms. The class centres can increase as well as decrease during the learning process. Moreover this algorithm easily fits in an interactive knowledge discovery application. The specification of GNG algorithm is given below:

Algorithm Growing Neural Gas; Parameters; age(a1,a2) : age of connection between two classes a1 and a2; agemax : maximal age of connection; εa : error of class a; εb ,εb : coefficient of adaptation of winning class and its neighbours; Procedure; 8.; 1.Initialize two classes A = {c1, c2}, t=0. Initialize the connection set.; 2.Select at random an input signal ξ.; 3.Determine the winner s1 and the second-nearest cluster s2 ,the closest to ξ.; 4.If a connection between s1 and s2 does not exist already, create it. Set the age of the connection to 0; C = C U {(s1,s2 )}.age(s1,s2) = 0; 5.Add the squared distance between the input signal and the winner to a local error variable: Δεs1 = II ξ - wsiII^2.; 6.Adapt the reference vectors of the winner and its direct topological neighbours by fractions : Δwsi = εb*(ξ -wsi) , Δwi = εn*(ξ -wn); 7.Increment the age of all edges emanating from si; 8.Remove edges with an age larger than amax. If this unit has no more emanating edges, remove the unit as well.; 9.If the number of input signals generated so far is an integer multiple of a parameter l, add a new unit r to the network and interpolate its reference vector from q and f, decrease the error variables of q and f.; 10.If a stopping criterion (e.g., net size or some performance measure) is not yet fulfilled, continue, return to step 2.

The third algorithm, called CURE (Clustering Using REpresentatives), is an agglomerative algorithm where disjoint clusters are successively merged until the number of clusters reduces to the desired number of clusters. CURE can identify clusters that are not spherical as well as clusters with wide variances in cluster size. These features are particularly interesting while clustering medical images. The specification of CURE algorithm is given below:

Algorithm CURE; Parameters; S: n voxels in d-dimensional space; k: number of clusters; u,v w, and x: clusters; u.rep: the set of representative points in u cluster; u.closest: the closest cluster to u; kd_tree: data structure that stores the representative points for every cluster; the heap: data structure that stores the entries for various clusters u arranged in the increasing order of the distances; Procedure cluster (S,k); 11.; 1.T:= build_kd_tree(S) /*all voxels are inserted into the k-d tree */; 2.Q := build_heap(S)/*each input voxel is considered as a separate cluster */; 3.while size(Q) "/> k do{; 4. u := extract_min(Q) /*extract the top element in Q */; 5. v := u.closest; 6. delete(Q,u); 7. w := merge(u,v) /*merge the closest pair of clusters u and v and compute new; representative point for the new merged cluster w which is inserted into T */; 8. delete_rep(T,u); delete_rep(T,v); insert_rep(T,w); 9. w.closest := x /* x is an arbitrary cluster in Q*/; 10. for each x Q do {; 11. if dist(w,x) < dist(w, w.closest); 12. w.closest := x; 13. if x.closest is either u or v {; 14. if dist(x,x.closest)< dist(x,w); 15. x.closest := closest_cluster (T,x,dist(x,w)); 16. else; 17. x.closest :=w; 18. relocate(Q,x); 19. }; 20. else if dist(x,x.closest) "/> dist(x,w) {; 21. x.closest := w; 22. relocate(Q,x); 23. }; 24. }; 25. Insert(q,w); 26. }

The algorithm CURE, in contrast to the previous algorithms, does not favor clusters with spherical shapes and similar volumes. The computational complexity of CURE is quadratic, so for large fMRI data it was necessary to employ random sampling of voxels, sacrificing clustering quality. More detailed discussion of CURE performance can be found in [8]

While the current approach supposes that the human expert builds up the hypothesis and the software, e.g. SPM [3], is only used to validate that hypothesis, data mining techniques can complement that approach by guiding the expert in his generation of new hypotheses, in particular by automatically showing up activated areas, and highlighting dependencies between those areas.

We have extended the SLICER system for clustering and interactive exploration of fMRI images features. In its current version, our software environment allows the physician to interact graphically with the clustering process, e.g. by modifying parameters of the algorithm, by focusing on a specific region of the brain, etc.

4. Visual data mining

The fMRI mining system integrates two functional parts: interactive clustering algorithms and visualisation package 3DSlicer developed by Harvard Medical School and AI Lab of the Massachusetts Institute of Technology (http://www.slicer.org [33]). The visualization package allows observing in 2D and 3D of the evolution of clusters discovered by data mining algorithms. The interface for interactive clustering has been also designed to engage the expert in the process of discovery. Therefore all clustering algorithms implemented in our system can be run in an interactive mode. The expert can follow the classification evolution from the beginning at regular time intervals.

Four levels of interaction can be distinguished. The lowest level corresponds to measurements provided on the state of classification and makes it possible for the user to make a decision when to interrupt the process. This interruption can be thus called upon systematically with regular intervals of times to take absolute measures on the classification and to visualize the intermediate results. This makes it possible to store or modify the current state of clustering. Thus, the expert may access the clustering algorithm and modify it at any point during its execution. Not only the parameters of the algorithm, but also the data space can be modified and refined during the clustering process. The second level concerns the definition and resizing of the data volume. The third level of interaction concerns the algorithms with the parameters that can be dynamically changed. Beyond these three levels of interaction, the management of saving and restoring the states of clustering constitutes the last level of classification interactions. Figure 2 illustrates the concepts of visual mining of fMRI data.

Three types of information are provided: statistical, temporal, and spatial. On the left part, the statistics describes the evolution of clusters; the user can examine the evolution of the errors, intra-class and inter-class inertias, the number of voxels connected to the center of each class, and the number of voxels which change class per unit time. The system provides also information that describes the dispersion of the clusters, the convergence of algorithm and the stability of clusters. Temporal information associated with each cluster represents graphically the fMRI signals over time. The signals are visualized in the form of curves on a graph. A paradigm, as in a sequence of stimuli and the model of the response, if available, can also be shown and compared online. On the right part, spatial information is simultaneously visualized in two manners: as 3D image, and as three perpendicular cuts across the brain. It is important to notice that during the clustering, the generated clusters are displayed in 2D and 3D. To facilitate interpretation, the 3D image of clusters can be superimposed, depending on the case, with a structural MRI of the subject or a standardized volume acting as an anatomical atlas.

Figure 2.

Interface of fMRI mining system

The user can easily optimize the clustering process by the online observation of statistical measures and their variations induced by modifications of parameters. Thus, the adjustment of the algorithm can be suitably carried out without continuously stopping and restarting the process.

To increase the speed of execution and to reduce the complexity of the results, the expert may find it beneficial to reduce the volume of the data. The space of exploration can be adjusted in several ways. The standard version of 3DSlicer assures not only the 3D visualization, but also allows the selection of volumes of interest or disinterest. The novel functionalities of the system permit the user also to restrict the space of research using anatomical structures or a subset of generated clusters. For instance, a threshold may eliminate the space that surrounds the brain. In the case of standardized images, it is possible to limit research to the grey matter. Thus, the expert can continue a classification by focusing his interest on certain regions after having eliminated voxels belonging to other non-relevant clusters.

5. Case studies

During the process of knowledge discovery, clusters generated by unsupervised classifiers must be validated by an expert, who retains only the relevant ones. Each selected cluster corresponds to a set of voxels, or zones of the brain, with a similar hemodynamic response over time. These responses can be explicitly characterized via the construction of classification rules. These rules combine observed temporal patterns with spatial information, such as the activity of voxels in neighbouring zones, or domain knowledge, such as the atlas of functional zones regions of the brain. It is important to note that the data normalisation (i.e. when the images are recalibrated to correspond to one brain type; voxels of several series of images of one person, or of several people), the voxels on the same positions correspond to the same zones of the brain. Temporal patterns can be synchronized with the paradigm, for example, to discover the interactions between areas regions of the brain used for visual memory. Temporal patterns can also be independent of any paradigm, for example, to highlight the succession of typical activations of region areas of the brain associated with hallucinations.

In the next part of the section, the performance of the clustering algorithms on the synthetic and real fMRI data will be described.

The synthetic data test

The synthetic data were composed of two parts, purely artificial simulated activations and real data. In the first experiment, the selected images corresponded to the auditive test conditions: "silence" and "talk". All the images of the "silence" condition were real data. Added to this series of 40 images were the synthetic activations formed by time series in crenel (square signal), simulating a paradigm of the block type. The localization of activations was a cubic volume of 5 voxels of each dimension. The level of average noise of these 125 voxels of the bottom was measured by taking double the variance of the intensities of the voxels over time. From this measurement, it was then possible to control the signal-to-noise ratio of the activations, by adding crenels of desired intensities to the considered voxels.

The performance evaluation was based on a confusion matrix indicating the true-positive ratios in particular. To create this matrix, it was necessary to determine the positive and negative classes. The figure 3 shows, in a simplified way, how the positive classes were defined and their coverage of a zone of activation.

Figure 3.

Simplified iIlustration of relation between a zone of activation and positive/negative classes

The table 1 illustrates the results of five clustering algorithms on synthetic data varying the parameter Signal-to-Noise (S/N). More detailed results may be found in research reports [Hommet, 2003].

S/N	GNG	SOM	LBG	K-means	CURE
1.2	55	0	0	15	52
1.4	100	0	15	35	76
1.6	100	70	55	44	72
1.8	100	100	65	50	72
2.0	100	100	69	70	72
3.0	100	100	100	82	88

Table 1.

Detection frequency (%) in respect of Signal-to-Noise ratio

Amongst implemented algorithms, two of them generated relatively stable and coherent clusters. The first, the Growing Neuronal Gas algorithm [7], has been adapted to the interactive classification of fMRI images. Contrary to the majority of the other methods, the originality of the GNG lies in the fact that the number of classes is not fixed in advance. In addition, the network topology and the number of neurons can be dynamically increased or decreased during the classification, making the algorithm efficient to cluster a large volume of data. The connections between neurons have the property of aging, and disappear when they reach a preset maximum age. This property is the cause of the disappearance of neurons, which are eliminated when they are not connected any more. The creation of the new neurons is made at regular time intervals by inserting a new neuron near the neuron with the greatest error. The error of a neuron is evaluated using the sum of the distances from this neuron to the voxels that declared this neuron as a winner. Thus the network is reinforced near the neurons which traverse the longest distances or which bring together the most voxels. Our contributions pertain to two aspects of classification, in particular reducing the volume of data and setting the parameters of the algorithm with the assistance of statistical tools and visual data mining.

Figure 4.

Results of clustering

The second algorithm, CURE [8], that has been also adapted to an interactive mode, generates clusters that have non-spherical forms with large variance, and, moreover, it allows a reduction in the processing time by aggregation and sampling fMRI data.

That the K-means had the worst level of detection is not surprising. K-means and LBG to some extent are strongly dependent on their initializations. But these two algorithms were not stripped of usefulness insofar as they did not require any parameter apart from the number of classes. This simplicity often enables them to obtain better results than GNG and SOM when the parameter setting of the latter is not optimal.

Block type data of a fMRI series of auditive tests

The real data used come from the site of the London research institute, "The Wellcome Department of Imaging neuroscience [http://www.fil.ion.ucl.ac.uk/spm/]; they form part of the test set of the program SPM99. The series was composed of 96 MRI acquisitions of the brain recorded with a repetition time of 7 s. The paradigm of the block type alternates the two following situations: a condition without stimulus and a condition of auditive stimuli consisting of repetitions of two-syllable words. With the pre-processed series, the classifications were performed by the four algorithms under the same conditions. With the rough series, one classification has been obtained using the paradigm data. The resulting classifications highlighted the most significant noises, such as the movements of the subject’s head, which generated clusters along the most intense contrasts of the image. By gathering the data in agreement with the paradigm, the zone activated by the auditive test could be revealed. The paradigm of this test was of the block type, where two conditions followed one another and were repeated, forming a periodic pattern. Here, the signals tend towards two conditions during each of the 6 images and repeat 8 times. Therefore, the 8 periods have been compressed into one period made up of 12 images. Classification of this small series preferentially reveals the expected zones in the form of compact blocks of voxels. However, within a class there are also scattered voxels. The treated data was compacted by a factor 16 compared to the pre-processed series. Amongst applied algorithms, two generated relatively stable and coherent clusters: GNG and SOM. Figure 4 illustrates clustering results of 4 independent runs of four algorithms (white coloured zones on fig.4 are relevant).

This result is interesting insofar as there was no pre-processing applied and the volume of the data was extremely reduced, considerably increasing the speed of execution.

6. Conclusion

In this article, a novel approach to interactive mining of fMRI data has been presented. The engagement of a physician in the process of knowledge discovery has been discussed, as well as the specificities of fMRI clustering with weak prior knowledge. Several clustering algorithms were evaluated. The experiments have shown that the Growing Neural Gas algorithm demonstrates the highest clustering performance and acceptable robustness.

The results have also demonstrated that the proposed new approach could be applied in detecting event-related fMRI data. It is important to underline that the parameters of the exploration algorithms can be modified during the course of execution. This dynamic aspect was a determining factor in its usage in interactive data mining.

The integration of clustering algorithms with the 3D Slicer has allowed visual exploration of created clusters and has provided the physician with more comprehensible information in quasi real time. The first results on synthetic data and block type data are encouraging and allow us to extend this work towards experiments with event-driven data where the signal-to-noise ratio is particularly weak and noisy data mask the relevant information.