Estimating the intrinsic dimension in fMRI space via dataset fractal analysis - Counting the `cpu cores of the human brain

Estimating the in trinsic dimension in fMRI space via dataset fractal analysis Coun ting the ‘cpu cores’ of the human brain Harris V. Georgiou (MSc, PhD) Dept. o f Informatics & T elecomm unications, National Kap o distrian Univ. of Athens (NKUA/UoA) Email: harris@xgeorgio.info –- URL: http://xgeorgio.info Abstract. F unctional Magnetic Resonance Imaging (fMRI) is a pow- erful non-in v asive tool for lo calizing and analyzing brain activity . This study fo cuses on one v ery imp ortant aspect of the functional prop erties of h uman brain, sp ecifically the estimation of the lev el of p ar al lelism when p erforming complex cognitive tasks. Using fMRI as the main mo dalit y , the human brain activit y is inv estigated through a purely data-driven signal pro cessing and dimensionalit y analysis approach. Specifically , the fMRI signal is treated as a multi-dimensional data space and its intrin- sic ‘complexity’ is studied via dataset fr actal analysis and blind-sourc e sep aration (BSS) metho ds. One simulated and tw o real fMRI datasets are used in com bination with Indep endent Comp onent Analysis (ICA) and fractal analysis for estimating the in trinsic (true) dimensionality , in order to provide data-driv en experimental evidence on the n umber of indep enden t brain pro cesses that run in parallel when visual or visuo- motor tasks are p erformed. Although this num b er is can not b e defined as a strict threshold but rather as a con tinuous range, when a sp ecific activ ation level is defined, a corresp onding num b er of parallel pro cesses or the casual equiv alent of ‘cpu cores’ can b e detected in normal human brain activity . Keyw ords: fMRI, ICA, fractal dimension, fractal analysis, human brain 1 In tro duction Human brain is the most adv anced and efficient signal-processing machine kno wn to da y . It corresp onds to only 2% of total b o dy weigh t in adults (ab out 1.5 kg), y et it consumes 20% of blo o d oxygen and 25% of glucose, with only 20W at p o wer p eak. It consists of roughly 100 billion neurons with 1,000-10,000 synapse in terconnections each, pack ed in 1130-1260 cm 3 of volume, making it the most complex organ in the human b o dy [42,14,30]. Analyzing its structure and func- tionalit y , esp ecially during the actual pro cess of some cognitive task or in relation to some mental impairment, has b een a scientific c hallenge for centuries. How- ev er, only recent technological breakthroughs hav e enabled the study of the inner 1. INTR ODUCTION w orkings of living brains. Ev en today , simulating the structure and only basic neuron functionalit y of a full-scale human brain in a digital computer is still an infeasible task. F unctional Magnetic Resonance Imaging (fMRI) [33,30,45] is a p ow erful non- in v asive tool for lo calizing and analyzing brain activity . Most commonly it is based on blo o d o xygenation level-dependent (BOLD) contrast, which translates to detecting lo calized changes in the hemodynamic flow of oxygenated blo o d in activ ated brain areas. This is ac hieved by exploiting the different magnetic prop erties of o xygen-saturated versus oxygen-desaturated hemoglobin. In the human brain, tasks in v olving action, perception, cognition, etc., are p erformed via the sim ultaneous activ ation of a n umber of functional br ain net- works (FBN), which are engaged in proper in teractions in order to effectively execute the task. Such netw orks are usually related to low-lev el brain functions and they are defined as a num b er of se gr e gate d sp ecialized small brain regions, p oten tially distributed o ver the entire brain. In order to properly detect these ac- tiv ations and iden tify the set regions that constitutes a FBN related to a sp ecific task, the 3-D space o ccupied by the brain is partitioned in to a grid of ‘cub es’ or voxels . Each vo xel constitutes the elementary spatial unit that acts as a signal generator, recorded and registered as a lo w-resolution 1-D time series. Actual fMRI vo xel signals from brain scans can b e considered as a mixture of v arious comp onen ts or sour c es with differen t temp oral and spatial characteristics. These sources can b e classified as of in terest and as artifacts [7]. In order to understand the true functionality and full p oten tial of the h uman brain, data-intensiv e approaches are required for analyzing the actual brain sig- nal (e.g. fMRI, EEG) during sp ecific cognitive tasks. Current research inv olves m ulti-disciplinary endeav ors, from Bio chemistry and Neurophysiology to Simu- lation and VLSI design, with pro jects like the Human Brain Pro ject (HBP) by EU [18] and Brain Researc h through A dv ancing Innov ative Neurotechnologies (BRAIN) by USA [47]. There is also v ery activ e researc h and developmen t effort in the industry , where pro jects like the recently announced ‘T rueNorth’ chip b y IBM implement a million-scale neural net work grid in sp ecial-purp ose VLSI and extremely high p o wer efficiency [37,44]. How ever, all these efforts are currently fo cused on the structur al prop erties of the human brain, i.e., the neural netw orks top ology and connectivit y , while the functional and higher-level cognitive prop- erties are still v ery difficult to mo del. In practice, this means that the hardware necessary to build and fully simulate (at the neuron cell level) an actual artificial ‘brain’ equiv alen t to a small animal’s now becomes av ailable, but the problem of turning this construction to a machine with actual cognitive and abstract func- tionalit y still remains (for the most part) unsolved, with only application-sp ecific mo dules b eing developed successfully (e.g. artificial retina implants, with some visual pro cessing capabilities [15]). This study fo cuses on one very imp ortan t asp ect of the functional prop erties of h uman brain, sp ecifically the estimation of the level of p ar al lelism when p er- forming complex cognitive tasks. In some very abstract sense, this is not muc h differen t than trying to recov er the (minim um) n um b er of actual ‘cpu cores’ 2 2. PR OBLEM DEFINITION required to ‘run’ all the active cognitive tasks that are registered in the entire 3-D brain volume while p erforming a typical fMRI exp erimental proto col that includes visual-only or visuo-motor tasks. Using fMRI as the main mo dalit y , the h uman brain activity is inv estigated through a purely data-driven signal pro cessing and dimensionality analysis ap- proac h. Sp ecifically , the fMRI signal is treated as a multi-dimensional data space and its in trinsic complexity is studied via dataset fr actal analysis and blind- sour c e sep ar ation metho ds. Section 2.1 provides an ov erview of the fMRI exp eri- men ts and the nature of sensory data; section 2.2 defines a prop er mathematical form ulation for the data unmixing task and its imp ortance in understanding the true sources of brain activity; section 3.1 provides hints to prop er data dimen- sionalit y reduction in fMRI; section 3.2 briefly describ es the basic metho dology for dataset fractal analysis and how it is applied for the estimation of the intrin- sic dimensionality of the fMRI data space; section 3.3 briefly describ es ICA as a t ypical approach for blind-source separation in signal pro cessing; sections 4.1, 4.2 and 4.2 describ e the simulated and real fMRI datasets used in this study; section 5 includes the exp eriments and results, using all the metho ds and datasets de- scrib ed earlier; sections 6 and 7 conclude the study with discussion of the results and their practical meaning. 2 Problem Definition 2.1 The nature of fMRI data In exp erimental fMRI pro cedures, there are tw o common activ ation schemes: the blo ck paradigms and the event-r elate d paradigms [5]. In the blo c k paradigm, the sub ject is presen ted with a sp ecific stimulus for a sp ecific time frame, e.g., a set of images of different placement, colors, patterns or categories, and the sub ject has to press a switch to signal positive or negative feedback as a resp onse. In the even t-related paradigm, the sub ject is exp osed to a series of randomized short-time inputs, e.g., a noise or a pain stimulus, with or without the need for sp ecific resp onse from the sub ject. In b oth cases, the external input is consid- ered as a primary ‘source’ and is temp orally correlated with the brain activity . Areas of high activ ation and correlation to the stimulation/response pattern are considered as highly relev ant to the sp ecific functional task (visual/motor cen- ters, pain receptors, etc). Th e same pro cedure can b e follow ed when there is no sp ecific external paradigm, constituting the ste ady-state functional analysis of brain activity . In this setup, there is no correlation to previously kno wn ac- tiv ation pattern and hence the analysis is essen tially a search for functionally indep enden t sources in the recorded fMRI signal. The acquired fMRI signal is registered in b oth spatial (3-D) and temp oral (1- D) domain, resulting in a comp osite 4-D signal. Each spatial axis is registered as a grid of spatial resolution 3-5 mm 3 , resulting in a 3-D grid of vo xels. T ypically , a complete v olume of vo xel data, e.g. 60x60x30 to 64x64x48, is recorded ev ery 1-2 seconds for a sequence of 100-150 time p oints [33,30,45]. This pro duces a total of roughly 108K-197K v oxels for ev ery time frame or, equiv alently , 11e6 3 2. PR OBLEM DEFINITION to 30e6 data p oints organized as a t w o-dimensional matrix, where eac h ro w corresp onds to a complete brain ‘snapshot’. In practice, the n umber of actual brain vo xels is smaller, since non-brain areas of the grid are masked out b efore an y further pro cessing; how ever, the data volume still remains within the same order of magnitude. A dditionally , typical fMRI exp erimental proto cols inv olve sev eral sub jects, in order to exclude any sub ject-sp ecific characteristics that may affect the statistical prop erties of the fMRI data under consideration. Clearly , this creates a high-volume data analysis pro cess that mak es it a very complex and computationally demanding task. In terms of signal pro cessing, the hemo dynamic resp onse function (HRF) [33,30,45] of the activ ated neurons, i.e., the changes in oxygen-ric h blo o d flow in the time domain, acts as a lo w-pass filter in the temp oral domain, which in turn mo difies the true activ ation signal that it is registered as fMRI data. In other w ords, the HRF of the activ ated neurons, i.e., the changes in oxygen-ric h blo o d flo w in the time domain, mo difies the true activ ation time-series signal that it is registered for eac h brain vo xel as fMRI data. Moreo ver, the HRF is kno wn to b e spatially-v arying, which means that there are slightly different hemo dynamic resp onses for differen t areas of the brain, as w ell as differen t HRF s betw een differen t sub jects. This means that traditional regression approac hes like General Linear Mo del (GLM) appro ximations [33,45,30] that require a pre-defined ‘design matrix’ are clearly sub-optimal, since it is typically constructed as p erm utations, transformations, time-shifts and deriv atives of one (assumed) ‘univ ersal’ HRF. There are also additional features that mak es this approximation even more difficult in practice, suc h as the fact that the vo xels’ activ ations are assumed to b e statistically indep endent (while lo cally they are not, due to the physical prop erties of the v eins and hemo dynamics), as well the v arious artifacts that are in tro duced to the signal by external factors (scanner drift, electronic noise, head mo vemen ts, respiration, cardiac pulsation, etc) [33,30,45]. 2.2 Understanding brain activity In fMRI analysis, the sources of interest include task-related, transiently task- related and function-related sources, meaning that in a task-sp ecific fMRI ex- p erimen t most of the task-related activit y is expected to be spatially isolated and temp orally synchronized with the corresp onding input/stimulation patterns. Therefore, these sources are exp ected to app ear as sup er-Gaussian in nature due to the spatial and temp oral lo calization of such task-related brain functionality . Sp ecial matrix factorization algorithms are required to reformulate the fMRI data as a multiplication of tw o other matrices, where one is for the time courses of the estimated signal ‘sources’ and one for the corresp onding spatial maps of related brain activity . F ormally put, if Y ∈ R t × n is the full fMRI data matrix with t rows as time p oin ts and n brain vo xels ‘unwrapped’ into a linear v ector, then the fMRI data matrix can be factorized as Y = T S , T ∈ R t × p , S ∈ R p × n , where the p spatial maps are collected as ro ws in S and eac h column of T con tains the activ ation pattern along time for the corresp onding spatial map. 4 3. DIMENSIONALITY ANAL YSIS OF THE FMRI DA T A SP A CE In GLM [33,45,30], T is the pre-defined ‘design matrix’ that contains p er- m utations, transformations, time-shifts and deriv atives of one (assumed) ‘uni- v ersal’ HRF, while S is the matrix of the corresponding regressor co efficients in eac h ro w. Although the GLM approach is sufficient w hen only sp ecific sensory- related signal sources (external stim uli) are considered, in the general case it is not p ossible to define a global design matrix for all signal sources and all (mul- tiple) sub jects. Instead, Indep endent Comp onent Analysis (ICA) is the most commonly used metho d for this task, in the con text of blind source separation (BSS) [25,8,27] (see section 3.3). While GLM and ICA are the dominating approaches for directed or blind unmixing mo dels, resp ectively , in fMRI analysis, the large v olume of v oxel data and the inherent properties of the fMRI signals make the unmixing task highly demanding in b oth memory and computational resources. Moreov er, prop er iden- tification of ‘universal’ FBNs requires multiple exp eriments with different sub- jects, which means working with m ultiplied volume of sensory data or combining m ultiple unmixing results ov er v arious runs [17,10]. In either case, unmixing al- gorithms are required to b e both fast and accurate in identifying the signal ‘sources’ of fMRI data and the activ ated areas in the brain corresp onding to the sp ecific paradigm source. Based on these prop erties, it is clear that the analysis of fMRI data, in the sense of its decomp osition in to distinct sources and the identification of the ones related to a sp ecific task or functional activity in the brain, is a v ery difficult task. The lac k of strict sp ecifications for a ‘univ ersal’ HRF and bac kground artifacts, hence in turn for an accurate pre-defined ‘design matrix’ for a standard GLM mo del, makes it a t ypical candidate for BSS approaches such as ICA, as well as alternativ e approac hes like Dictionary Learning (DL) and Compressive Sensing (CS). Recently , there is an increased interest for alternativ es to ICA for data- driv en fMRI unmixing. Notably go o d results hav e b een attained with Dictionary Learning (DL) - based fMRI analysis [32,2,3]. Also, an improv ed v ariation of K- SVD was prop osed as the basis for Dictionary Learning (DL), customized to fMRI analysis [29]. 3 Dimensionalit y an alysis of the fMRI data space 3.1 Data decimation and intrinsic dimensionalit y One wa y to deal with the high complexity of the BSS task in fMRI data is to reduce the num b er of vo xels under consideration. Sp ecifically , adjacen t neurons in the brain can b e considered highly correlated in terms of their resp onses to external stimuli, provided that the blo o d vessel netw orks at very small scales ac- tually introduce some spatio-temp oral correlation. Hence, their BOLD resp onse and HRF can b e considered, at some degree, statistically de pendent. If the spa- tial resolution of the fMRI signal is high, adjacent v oxels in the original 2-D or 3-D volume scan can b e considered statistically redundant. Therefore, some form of decimated vo xels set can b e used instead as input for the unmixing task, 5 3. DIMENSIONALITY ANAL YSIS OF THE FMRI DA T A SP A CE without sacrificing the accuracy of identifying the true inherent sources of the data. In fMRI, using decimated v ersions of the original data in the spatial and/or the temporal domain is not an uncommon practice. Indeed, some w orks refer to sub-sampling the fMRI in a spatial (vo xels) or temp oral (time p oin ts) sense, when the resolution is considered high enough. How ever, this has b een applied only as a one-time pre-pro cessing step in the preparation of data, i.e., b efore any real BSS analysis is conducted (e.g. see in [31,32]). F urthermore, the decimation ratio used in eac h case is chosen in a purely empirical wa y , since there are cur- ren tly no analytical studies with regard to the resulting qualit y of the decimated fMRI data. Spatio-temp oral correlations b etw een v oxels and statistical dep endencies are essen tially the reason why the fMRI data space has an intrinsic (true) dimen- sionalit y muc h smaller than the num b er of v oxels, i.e., the data matrix Y ∈ R t × n is of column rank c  n . How ever, for prop er unmixing of the fMRI data, the column rank of matrix Y should b e retained even when some decimation pro cess is employ ed. In other words, the selection of a smaller subset of vo xels (instead of all) should b e conducted in a wa y that do es not destroy the information con tent of the full data, but instead exploit the the fact that the num b er of vo xels n is v ery large and their inherent statistical prop erties can b e prop erly retained with a m uch smaller subset. In the cases when only a small set of the signal sources are considered, i.e., the time series of some external stimuli (plus some transformations of it), then regression metho ds like GLM can b e easily formulated with the prop er ‘design matrix’ to recov er the related brain activit y . When the analysis is conducted in the BSS sense, i.e., all ma jor signal sources are to b e reco vered (including the stimuli time series), then decomp osition methods lik e ICA pro vide a w ell- form ulated statistical framework for this task, as long as the prop er constraints are asserted as v alid (most imp ortan tly , the assertion of at most one Gaussian signal source). How ever, when these statistical assertions are not fully satisfied or when there is a large num b er of signal sources that are ‘exp onentially decay- ing’ in terms of imp ortance (con tribution to the mixed signal’s v ariance, p ow er sp ectrum and approximation error), then the num b er of indep endent compo- nen ts that ICA or other similar algorithms is limited only b y some external pre-defined threshold. In other words, the data matrix Y ∈ R t × n can b e factor- ized appr oximately as Y ' T S , T ∈ R t × p , S ∈ R p × n , with the reconstruction error b ecoming smaller as p increases. In theory , if the true sources of the mixed signal are p erfectly separable in the BSS sense, then ICA will stop after recov- ering exactly p = c comp onents, where c  n is the column rank of the data matrix Y . This means that there are exactly p comp onents, i.e. time courses and corresp onding activ ation maps, that can fully reconstruct the fMRI data for the en tire brain activity . Hence, the definition of the optimal v alue for p b y means of non-parametric (data-driven) estimation pro cedure is of utmost importance in the BSS task for fMRI unmixing. 6 3. DIMENSIONALITY ANAL YSIS OF THE FMRI DA T A SP A CE 3.2 Dataset fractal analysis In recent years, dimensionality analysis in signal pro cessing has b een extensively link ed to fractal analysis and fr actal dimension , as a non-parametric metho d for the quantitativ e characterization of the complexity or ‘randomness’ of a signal [16,34]. When applied to 1-D signals, metrics like the Hurst exp onent or Lya- punov exp onent hav e b een used as statistical features to describ e v arious t yp es of data series, from biomedical signals (e.g. EEG, ECG, etc) to financial and climate time series. In 2-D signals, these metho ds provide additional features for characterizing the texture of images, e.g. when analyzing biomedical modal- ities (radiology , ultrasound, MRI, etc) [9]. F ractal dimension is closely linked to these fractal parameters and it provides a clear distinction b etw een the emb e d- ding space, i.e., the full-rank space in the algebraic sense, from the actual space spanned b y the registered sensory data. In the general case when fractal analysis is applied to some multi-dimensional signal, the estimation of the fractal dimen- sion can b e used as a realistic ev aluation of the ‘complexity’ of the space spanned b y the actual data points av ailable and, hence, a v ery useful hint regarding the inheren t redundancy in a given dataset. In order to establish a preliminary estimation of the complexit y and intrinsic dimensionalit y of datasets, fractal analysis provides a data-cen tric approach for this task. Dataset fractal analysis, sp ecifically the calculation of intrinsic fr actal dimension (FD) of a dataset, provides the quantitativ e means of inv estigating the non-linearit y and the correlation b etw een the av ailable ‘features’ (i.e., di- mensions) by means of dimensionalit y of the em b edding space [9,39]. F ractal dimension has also b een used as an alternative w ay of c haracterizing the dis- criminativ e pow er of eac h ‘feature’ separately , th us pro viding a non-statistical w ay of ranking them in terms of imp ortance, e.g. as means of non-parametric feature selection in classification tasks [46]. The fractal analysis of datasets has b een used successfully in previous studies [35,19,36] and it has b een prov en very v aluable as a to ol for comparing arbitrary datasets of extracted features with the qualitativ e clinical prop erties that an exp erience physician uses to c haracterize a mammographic image. The tw o most commonly used metho ds of calculating the fractal dimension of a dataset are the p air-c ount ( P C ) and the b ox-c ounting ( B C ) algorithms [34,39,43,16]. In the pair-count algorithm, all Euclidean distances betw een the samples of the dataset are calculated and a closure measure is then used to cluster the resulting distances space into groups, according to v arious ranges r , i.e., the maximum allow able distance within samples of the same group. The P C v alue is calculated for v arious sizes of r and it has b een pro ved that P C ( r ) can b e appro ximated by: P C ( r ) = K · r D (1) where K is a constant and D is called pair-count exp onen t. The P C ( r ) plot is then a plot of: l og ( P C ( r )) versus l og ( 1 / r ) , i.e., D is the slop e of the linear part of the P C ( r ) plot ov er a sp ecific range of distances r . The exponent D is called c orr elation fr actal dimension of the dataset, or D 2 . 7 3. DIMENSIONALITY ANAL YSIS OF THE FMRI DA T A SP A CE The box-coun ting approach calculates the exp onen t D in a sligh tly different w ay , in order to accommo date case of large datasets with size in the order of thousands; ho wev er, it essentially calculates an approximation of that same cor- relation fractal dimension v alue, i.e., D 2 . It is commonly used when the datasets con tain large num b er of samples, usually in the order of thousands [46,4]. In this case, instead of calculating all distances betw een the samples, the input space is partitioned into a grid of n -dimensional cells of side equal to r . Then, the samples inside each cell are calculated and the frequency of o ccurrence R r , i.e., the count of samples in a cell, divided by the total num b er of samples, is used to appro ximate the correlation fractal dimension by: D 2 = ∂ l og P i  R i r  2 ∂ l og ( 1 / r ) (2) Ideally , b oth pair-coun t algorithm and b o x-counting algorithm calculate the same v alue, i.e., the correlation fractal dimension D 2 of the initial dataset, which c haracterizes the in trinsic (true) dimension of the input space [4]. In other words, D 2 w ould b e the minimum dimension of the dataset if only ‘p erfect’ features w ere allow ed, i.e., totally uncorrelated and with the b est discriminative p o wer a v ailable within the sp ecific set of features. In this study , fractal feature analysis w as applied to b oth the initial set of qualitativ e characteristics, provided by the exp ert physician, as w ell as the con- structed datasets of morphological features. In all cases, the pair-count algorithm emplo ying Euclidean distances was used, due to the relativ ely small num b er of samples av ailable, as well as the b etter stability and accuracy for D 2 against the b o x-counting approach [43]. In order to calculate the slop e at the linear part of the P C ( r ) plot, a para- metric sigmoid function w as used for fitting betw een the sample p oin ts of the plot. In the parametric sigmoid function: y = y 0 + C y  1 1 + exp ( − C x ( x − x 0 ))  (3) the ( x 0 , y 0 ) identifies the transp osition of the axes, while C x and C y iden tify the appropriate scaling factors. Sp ecifically , the v alue of C x affects the steepness of the central part of the curve, while C y sp ecifies the Y -axis width of the sigmoid curv e. Then, the slop e of the linear part around the central curv ature p oint, i.e. the v alue of D 2 , is: ∂ 2 y ( x 0 ) ∂ x 2 = 0 ⇒ D 2 = ∂ y ( x 0 ) ∂ x = C x · C y 4 (4) The fitness of the parametric sigmoid ov er a range of samples assumes uni- form error weigh ting o ver the en tire range of data. Thus, if a large p ercen tage of p oin ts lies near the upp er bound ( y = y max ) or low er b ound ( y = y min ) of the Y -axis range, as in most cases of P C ( r ) plots, then the fitness in the central region of the sigmoid, i.e., where the slop e is calculated, can b e fairly p o or. F or 8 3. DIMENSIONALITY ANAL YSIS OF THE FMRI DA T A SP A CE this reason, an additional weigh ting factor was introduced in the fitness calcu- lation in this study . Sp ecifically , the T ukey (tap ered cosine) parametric window function [21] was applied ov er the Y -axis range when calculating the ov erall fit- ness error of the sigmoid. The T ukey window is parametric ( q -v alue) in terms of the exact form around its cen ter, ranging from completely rectangular ( q = 0 ) to completely triangular or Hanning windo w ( q = 1 ) . When applied o v er the Y -axis range, the rectangular case is equiv alent to calculating the fitness error uniformly ov er the entire range, while the triangular case is equiv alen t to calcu- lating the fitness error primarily against the central p oin t of the sigmoid curv e. In this study , all fitness calculations employ ed T ukey windows as error w eight- ing factors, using parameters q in the range b etw een 0.5 and 1.0 for optimal slop e results. The equation for computing the co efficien ts w j of a discrete T ukey windo w of length N ( j = 1 ...N ) is as follows: w j =        1 2  1 + cos  2 π ( j − 1) q ( N − 1) − π  , 1 ≤ j < q 2 ( N − 1) 1 , q 2 ( N − 1) ≤ j ≤ N − q 2 ( N − 1) 1 2  1 + cos  2 π q − 2 π ( j − 1) q ( N − 1) − π  , N − q 2 ( N − 1) < j ≤ N (5) 3.3 Indep enden t Comp onen t Analysis (ICA) In blind source separation (BSS), ICA has b een successfully applied to fMRI data for many years [25,8,27,13]. Since the fMRI consists of a mixture of unknown comp onen ts, corresp onding to different brain sources of activity , the unmixing pro cedure is essen tially a BSS problem. How ev er, due to the relatively low tem- p oral and spatial resolution of fMRI data, the non-stationary prop erties of the signal due to brain- and mac hine-state v ariations, as well as the unknown num- b er and exact statistical prop erties of the sources, the BSS of fMRI data is not a trivial task. ICA is based on iden tifying non-Gaussian prop erties betw een the sources and separating them from the mixture, essentially reconstructing the original signal as a linear com bination of iden tified components, i.e., similarly to the previously discussed formulation Y = T S , T ∈ R t × p , S ∈ R p × n . In this case, S is the matrix of indep enden t components (spatial maps of brain activity) and T is the mixture matrix (corresp onding time courses). In fMRI, the ICA can be p erformed in the spatial or temp oral dimension of the (v ectorized) vo xel data, pro ducing either spatial or temporal components in matrix S . Several studies ha ve b een conducted in whether spatial or temp oral ICA works b etter for BSS in fMRI data [8]; how ever spatial maps, i.e., retrieving S as spatial comp onen ts, seem to b e more accurate and useful in most clinical applications of fMRI. The t wo most common approaches for ICA are the Infomax [6] and fastICA [ 26,24,25] algorithms. Although ICA has b een widely studied and emplo yed in fMRI, recen t works ha ve identified the relev ant adv antages of analyzing brain activity under the 9 4. D A T ASETS sparsit y , instead of statistical indep endency , of the underlying mixture of indi- vidual comp onents. Additionally , the BSS problem itself has b een identified b y few researchers as equiv alent to Dictionary Learning (DL) [1,2,32,31,29], which is already used in v arious applications. Since ICA do es not include an y sparsit y constrain ts (lik e in DL), while at the same time it assumes sp ecific statistical prop erties for the underlying signal sources (at most one Gaussian distribution, minimal noise artifacts). Hence, ICA unmixing of fMRI data that do not fully satisfy these constraints will construct factorizations that include the maximum allow able n umber of comp onents for the reconstruction of the original (mixed) data with the minimum error. In other w ords, as describ ed in section 3.1, when the fMRI data include non-trivial mix- tures of sources (as in the case of the simulated dataset, see section 4.1), ICA will construct a factorization mo del Y ' T S , T ∈ R t × p , S ∈ R p × n , with p = p max and non-zero reconstruction error. Similar problems emerge when using sparsity- a ware approaches as in DL [29], since they t ypically pro duce factorizations with p  p max (here, p max is the dictionary size), but with larger reconstruction errors, as exp ected. In this study , ICA is used as one of the most p opular approaches in BSS problems like the fMRI unmixing task. In the sim ulated fMRI datasets (see sec- tion 4.1), ICA provides an exact estimation of the intrinsic dimensionality of the signal, which is exp ected to b e lo w er than the pre-defined sources used in the mixture. In the real fMRI datasets (see sections 4.2 and 4.2), ICA pro vides appro ximate factorization mo dels and a quan titative wa y to track the signal re- construction error as the num b er of used comp onen ts changes. In b oth cases, the factorization mo dels pro vided b y ICA are used as a v erification tool for v alidat- ing the quality of the estimated fractal dimension of eac h dataset. Although the exact num b ers b etw een p max and the fractal dimension calculated differ due to their inherently different meaning, trac king their changes in parallel and com- paring results is used here as a v aluable to ol for dimensionality analysis of the fMRI datasets. 4 Datasets The inv estigation of fMRI space complexity and in trinsic dimensionality was conducted with tw o separate types datasets, namely one of simulated fMRI data and tw o of real fMRI data from carefully designed exp eriments. The sim ulated data were introduced as the means to v erify the recov erability of the intrinsic dimension when the real signal sources are known and w ell-defined, while the real data were used as guidelines for estimating the true brain activity in tw o t ypical cognitive tasks (visual recognition task and visuo-motor task). 4.1 Sim ulated fMRI datasets In this study , an adapted version of the real-v alued fMRI data generator co de from the MSLP-Lab [38] to olb ox was used for creating artificial fMRI data as a 10 4. D A T ASETS mixture of eight main sources [13]. Using the basic kno wledge of the underlying statistical characteristics of the underlying sources, the comp onents include three highly sup er-Gaussian sources (S1, S2, S5), a Gaussian source (S4) and a sub- Gaussian source (S3), plus tw o more sup er-Gaussian sources (S6, S8) and a sub- Gaussian source (S7). The time course for eac h comp onent defines the temp oral c haracteristics of the corresp onding source, namely one task-related (S1), tw o transien tly task-related (S2, S6) and several artifact types (S3, S4, S5, S7, S8), including respiration, cardiac pulsation, scanner drift, bac kground noise, etc. These sources can b e considered as spatial maps that are activ ated according to their time course and mixed linearly to pro duce the final (simulated) fMRI data. Although in typical fMRI exp erimen ts there is only one sensory ‘input’ (stim- ulation), here the full set of eigh t sources (S1...S8) w as considered throughout the ev aluation. Sp ecifically , the sim ulated fMRI data included eight spatial maps of size 60x60 vo xels (2-D ‘slices’) and a 100-p oint time course, with statistical prop- erties as describ ed ab ov e. Each spatial map was linearized by row-concatenation in to a (row) vector of 3600 vo xels, registered along its time course (column) v ector of 100 p oints. Finally , these eight 100x3600 matrices of spatio-temp oral maps were mixed linearly to pro duce the final eight-source mixing of simulated fMRI data in to one matrix of that same size. Hence, in terms of the problem form ulation presented in section 2.2, the final matrix of (simulated) fMRI data is registered as Y ∈ R t × n , where t = 100 time p oints and n = 60 2 = 3600 v oxels. Since the final data matrices are alwa ys linearized in a similar wa y b efore applying any unmixing algorithm, using 2-D ‘slices’ of (simulated) vo xels instead of full 3-D (real) brain scans in each time p oint affects only the volume of the data and not the task itself. 4.2 Real fMRI datasets ds101 – The ‘Simon’ task The ‘NYU Simon T ask’ dataset [28] comprises of data collected from 21 healthy adults while they p erformed a rapid ev ent-related Simon task. On each trial, the inter-trial interv al (ITI) w as 2.5 seconds, with n ull ev ents for jitter), a red or green b ox app eared on the righ t or left side of the screen. Participan ts used their left index finger to resp ond to the presen tation of a green b o x, and their right index finger to resp ond to the presentation of a red b o x. In congruent trials the green b ox app eared on the left or the red b o x on the righ t, while in more demanding incongruen t trials the green b o x app eared on the righ t and the red on the left. Sub jects performed tw o blo c ks, eac h con taining 48 congruent and 48 incongruen t trials, presented in a pre-determined order (as p er OptSeq), in tersp ersed with 24 null trials (fixation only). F unctional imaging data were acquired using a research dedicated Siemens Allegra 3.0 T scanner, with a standard Siemens head coil, lo cated at the NYU Cen ter for Brain Imaging. The data obtained were 151 contiguous echo planar imaging (EPI) whole-brain functional volumes (TR=2000 ms; TE=30 ms; flip angle=80, 40 slices, matrix=64x64; FO V=192 mm; acquisition vo xel size=3x3x4 mm 3 ) during each of the tw o Simon task blo c ks. A high-resolution T1-weigh ted anatomical image was also acquired using a magnetization prepared gradient 11 4. D A T ASETS ec ho sequence (MPRAGE, TR=2500 ms; TE= 3.93 ms; TI=900 ms; flip angle=8; 176 slices, F OV=256 mm). The ‘ds101’ dataset is hosted by Op enfmri.org for public access in ra w NiFTI format [40], including vo xel masks and brain map template, but without any pre- pro cessing (head mo vemen t, sensor drift, etc). In this study , the data from nine (out of 21) sub jects were used, including t wo runs each, for a total of 18 datasets of fMRI scans. Each dataset w as masked for exclusion of non-brain areas and thresholded for exclusion of brain areas with near-zero activity . The resulting num b er of vo xels ranged roughly b et ween 28K and 39K, while the num b er of snapshots was fixed to 151 time p oin ts. In terms of the form ulation of section 2.2, each fMRI data matrix is Y ∈ R t × n with t = 151 time p oin ts and 27631 ≤ n ≤ 38735 ‘non-zero’ vo xels. Three v ariants of each dataset w ere used, regarding the smo othing pre- filtering. Sp ecifically , according to standard fMRI acquisition practice, a Gaus- sian smo othing kernel was applied to the original 3-D vo xel space, in order to sup- press noise artifacts and enhance the spatial contin uity of the vo xel data. With resp ect to their F ul l Width at Half Maximum (FWHM) [11,12], or 2 √ 2 · l n 2 · σ ' 2 . 35482 · σ for Gaussian kernels, tw o different spatial sizes were used: 4 mm 3 and 8 mm 3 . In practice, since the vo xel resolution in this dataset is 3x3x4 mm 3 , the smaller k ernel p erforms (softer) av eraging on 1-1.33 neighboring vo xels, while the larger kernel p erforms (more aggressive) av eraging on 2-2.67 neighboring vo xels. These tw o ‘smo othed’ versions, plus the original non-smo othed version, are the three v ariants of each dataset, used throughout the exp eriments (see section 5 for details). ds105 – Visual ob ject recognition task The ‘Visual Ob ject Recognition T ask’ dataset [22,20,23,41] comprises of data collected from six healthy adults while they p erformed a visual recognition task. Neural resp onses, as reflected in hemo dynamic c hanges, were measured in six sub jects (five female and one male) with gradien t ec ho echoplanar imaging (EPI) on a GE 3T scanner (General Electric, Milwauk ee, WI) (rep etition time (TR) = 2500 ms, 40 3.5-mm-thic k sagittal images, field of view (FO V) = 24 cm, echo time (TE) = 30 ms, flip angle = 90), while they p erformed a one-back rep etition detection task. High- resolution T1-weigh ted sp oiled gradient recall (SPGR) images were obtained for eac h sub ject to pro vide detailed anatom y (124 1.2-mm-thick sagittal images, F OV = 24 cm). Stim uli were gray-scale images of faces, houses, cats, b ottles, scissors, sho es, c hairs, and nonsense patterns. The categories were chosen so that all stim uli from a given category would hav e the same base level name. The sp ecific categories w ere selected to allow comparison with our previous studies (faces, houses, chairs, animals, and to ols) or ongoing studies (sho es and b ottles). Control nonsense patterns were phase-scrambled images of the intact ob jects. T welv e time series w ere obtained in each su b ject. Each time series b egan and ended with 12 s of rest and contained eight stimulus blo c ks of 24-s duration, one for eac h category , separated by 12-s in terv als of rest. Stimuli w ere presen ted for 500 ms with an 12 5. EXPERIMENTS AND RESUL TS in terstimulus interv al of 1500 ms. Rep etitions of meaningful stimuli w ere pictures of the same face or ob ject photographed from different angles. Stim uli for each meaningful category w ere four images each of 12 different exemplars. The ‘ds105’ dataset is hosted by Op enfmri.org for public access in ra w NiFTI format [40], including vo xel masks and brain map template, but without any pre- pro cessing (head mo vemen t, sensor drift, etc). In this study , the data from six (all) sub jects were used, including three (out of 12) runs each, for a total of 18 datasets of fMRI scans. Each dataset was mask ed for exclusion of non-brain areas and thresholded for exclusion of brain areas with near-zero activit y . The resulting num b er of vo xels ranged roughly b et ween 22K and 47K, while the num b er of snapshots w as fixed to 121 time p oin ts. In terms of the form ulation of section 2.2, each fMRI data matrix is Y ∈ R t × n with t = 121 time p oin ts and 22387 ≤ n < 47192 ‘non-zero’ v oxels. Three v ariants of each dataset w ere used, regarding the smo othing pre- filtering. Sp ecifically , according to standard fMRI acquisition practice, a Gaus- sian smo othing kernel was applied to the original 3-D vo xel space, in order to sup- press noise artifacts and enhance the spatial contin uity of the vo xel data. With resp ect to their F ul l Width at Half Maximum (FWHM) [11,12], or 2 √ 2 · l n 2 · σ ' 2 . 35482 · σ for Gaussian kernels, tw o different spatial sizes were used: 4 mm 3 and 8 mm 3 . In practice, since the vo xel resolution in this dataset is 3x3x3.5 mm 3 , the smaller kernel p erforms (softer) av eraging on 1.14-1.33 neigh b oring vo xels, while the larger kernel p erforms (more aggressive) av eraging on 2.29-2.67 neighboring v oxels. These tw o ‘smo othed’ versions, plus the original non-smo othed v ersion, are the three v ariants of each dataset, used throughout the exp eriments (see section 5 for details). 5 Exp erimen ts and Results T wo separate sets of experiments were conducted in this study , one for BSS unmixing via ICA and one for dataset fractal dimension estimation. Both sets included all three fMRI datasets, namely one of simulated fMRI data and tw o of real fMRI data exp erimen ts (see sections 4.1, 4.2, 4.2). 5.1 ICA analysis of the datasets The ICA experiments that w ere conducted with the simulated fMRI data in- cluded tw o distinct realizations of the dataset, generated b y the same pro cedure and the same sp ecifications as describ ed in section 4.1. Since the data generation includes several noise components, the tw o realizations were used as an additional v erification chec k to v alidate that slightly different mixtures of (artificial) fMRI data do not pro duce significant differences in the ICA error-versus-components plots and estimated dataset fractal dimension. Figure 1 presents the time courses of the ICA factorization (matrix T ), with the blue curves representing each of the eight ideal (true) sources and the red curv es representing the corresp onding ICA-recov ered sources. Figure 2 illustrates 13 5. EXPERIMENTS AND RESUL TS the corresp onding activ ation maps (matrix S ) reco vered by ICA, spatially re- shap ed into prop er 2-D brain ‘slices’, where the reconstruction errors are visible as artifacts (‘ghost’ artifacts). Fig. 1. Ideal (blue) and ICA-recov ered (red) time courses of the eight sources in the sim ulated fMRI dataset. Parameter r is the correlation coefficient b etw een the original (ideal) and the recov ered time course, p is the corresponding p-v alue and rmse is the matc hing error. The first comp onent (upp er-left corner) corresp onds to the pre-defined external stimuli. Figure 3 presen ts the plot of reconstruction error (RMSE) v ersus the num b er of ICA comp onents used. Sp ecifically , after the ICA unmixing mo del is complete, the ICA comp onen ts are used one by one in rank-1 reconstructions of the orig- inal data and the corresp onding errors are used for sorting the comp onents in ascending order (smallest RMSE first). Subsequen tly , a set of comp onents starts from the first one (top of the list) and increased by adding the next one in each step, while estimating and registering the corresp onding reconstruction error. The plot illustrates the total reconstruction error decreasing almost line arly as the num b er of used comp onen ts increases, as exp ected. It should b e noted that for ‘p erfect’ ICA factorizations, as in the case of simulated fMRI data, the num- b er of comp onen ts recov ered by ICA is exactly the same as the num b er of signal sources (true) used in the mixture that created these data (see section 4.1). The ICA exp eriments that were conducted with the real fMRI data included t wo distinct datasets, ‘ds101’ and ‘ds105’, as describ ed in sections 4.2 and 4.2, resp ectiv ely . Instead of a single 2-D brain ‘slice’ as in the case of the simulated fMRI data, here the datasets employ full 4-D fMRI data, i.e., 3-D vo xel grid of the brain volume evolving in 1-D time course. Figure 4 illustrates a real example of a 2-D brain ‘slice’ for a single time p oin t, as it is registered in the ‘ds101’ dataset; the data are in raw unpro cessed mo de and no bac kground-exclusion masking. 14 5. EXPERIMENTS AND RESUL TS Fig. 2. ICA-reco vered activ ation maps of the eigh t sources in the simulated fMRI dataset, spatially reshaped into proper 2-D brain ‘slices’. The lo wer-left box corre- sp onds to the activ ation areas for the pre-defined external stimuli. The low er-right b ox illustrates the complete reconstructed fMRI mixture at time p oint t = 150 . Fig. 3. Reconstruction error versus num b er of used comp onents. ICA detects exactly eigh t comp onents, i.e., the num b er of true signal sources in the original mixture, and the final reconstruction error is practically zero. 15 5. EXPERIMENTS AND RESUL TS Signal-indep enden t noise is evident in the flat/blue areas, i.e., are outside the brain v olume. Fig. 4. Real example of a 2-D brain ‘slice’ for a single time p oin t, as it is registered in the ‘ds101’ dataset (unprocessed, no background-exclusion masking). Using the GIFT to olb ox for Matlab [13], Figure 5 illustrates the ICA-recov ered time course (red plot) and the corresp onding 2-D ‘flattened’ activ ation map that represen ts the actual response of the human brain in a visuo-motor task v ery similar to the exp erimental proto col employ ed in the ‘ds101’ dataset. Here, the ICA successfully reco vered one particular comp onent very similar to the exter- nal stimuli, which ideally is a square-shap ed pulse mo dulated b y the HRF (see section 2.1), instead of the noisy sin usoid curve. Figure 6 illustrates 10 of the 50 ICA-recov ered time courses of comp onents in a sample run with the ‘ds101’ dataset. Although the ICA conv erged successfully with the minimum attainable reconstruction error, the unmixing mo del failed to detect one single component that closely matches the ideal time course of the stim uli. Ho wev er, it is evident that one comp onen t (third from top-left) matches comp onen t no.7 and t wo comp onents (upp er/low er left) match comp onent no.8 of the simulated fMRI data as illustrated in Figure 1 in terms of ov erall shap e and noise c haracteristics. With respect to reconstruction error v ersus num b er of used components, Figure 7 and Figure 8 illustrate how the RMSE changes (drops) as the size of the ICA mixture b ecomes larger. Red curv es represen t the RMSE against the n umber of used comp onen ts up to an upp er limit of 10, 25, 50 and 100. The final 16 5. EXPERIMENTS AND RESUL TS Fig. 5. Sample result from the GIFT to olb o x for Matlab [13], illustrating the ICA- reco vered time course (red plot) and the corresp onding 2-D ‘flattened’ activ ation map of the actual resp onse of the human brain in a visuo-motor task very similar to the exp erimen tal proto col employ ed in the ‘ds101’ dataset. 17 5. EXPERIMENTS AND RESUL TS Fig. 6. ‘ds101’ (non-smo othed), 10 of the 50 ICA-reco vered time courses of comp onents in a sample run. (righ t-most) p oint in blue represents the maxim um-size, low est-RMSE in each case. Hence, the general slop e of the red curves, as well as the dotted blue line connecting the end p oin ts, illustrate the robustness of the ICA unmixing pro cess in eac h of the real fMRI dataset. 5.2 Dataset fractal analysis and intrinsic dimensionalit y Similarly to the ICA exp erimen ts, the dataset fractal analysis that was conducted with the sim ulated fMRI data included tw o distinct realizations of the dataset, generated by the same pro cedure and the same sp ecifications as describ ed in section 4.1. Since the data generation includes several noise comp onents, the t wo realizations were used as an additional verification chec k to v alidate that sligh tly different mixtures of (artificial) fMRI data do not produce significant differences in the estimation of the fractal dimension of the dataset. Figure 9 presents the log-log plot for the b ox-coun ting metho d of estimating the fractal dimension (FD) in the simulated fMRI dataset, as describ ed in sec- tion 3.2. Sp ecifically , the blue p oints represent instances of l og ( P C ( r )) v ersus log ( 1 / r ) , and the blue curv e is the b est-fit parametric sigmoid that is describ ed b y Eq.3. The FD is recov ered as the slop e of the curve in the central p oint at ( x 0 , y 0 ) , according to Eq.4. Based on the b ox-coun ting approach describ ed in section 3.2, the tw o real- izations of the simulated fMRI data resulted in FD v alues of 3.774 and 3.884, using the complete dataset with no decimation (100 sample vectors). This means that an av erage mean v alue of 3.83 can b e considered as a reliable estimate of the fractal dimension of the dataset. 18 5. EXPERIMENTS AND RESUL TS Fig. 7. ‘ds101’ (non-smo othed), ICA reconstruction error versus num b er of used com- p onen ts. Red curves represent the RMSE against the n umber of used comp onen ts up to an upp er limit of 10, 25, 50 and 100. The final (right-most) p oint in blue represents the maximum-size, low est-RMSE in each case. Fig. 8. ‘ds105’ (non-smo othed), ICA reconstruction error versus num b er of used com- p onen ts. Red curves represent the RMSE against the n umber of used comp onen ts up to an upp er limit of 10, 25, 50 and 100. The final (rightmost) p oint in blue represents the maximum-size, low est-RMSE in each case. 19 5. EXPERIMENTS AND RESUL TS Fig. 9. The log-log plot l og ( P C ( r )) v ersus log ( 1 / r ) (blue points) and the best-fit parametric sigmoid (red curv e) that reco vers the fractal dimension (‘FDE’) of the sim ulated fMRI dataset. The dataset fractal analysis exp eriments that were conducted with the real fMRI data included t wo distinct datasets, ‘ds101’ and ‘ds105’, as described in sections 4.2 and 4.2, resp ectively . Instead of a single 2-D brain ‘slice’ as in the case of the simulated fMRI data, here the datasets employ full 4-D fMRI data, i.e., 3-D v oxel grid of the brain volume ev olving in 1-D time course. Figure 10 presents the log-log plot for the b ox-coun ting metho d of estimating the fractal dimension (FD) in the ‘ds105’ dataset, as describ ed in section 4.2; the corresp onding plot for ‘dsd101’ is similar (omitted). Sp ecifically , the blue p oin ts represen t instances of l og ( P C ( r )) versus l og ( 1 / r ) , and the blue curve is the b est-fit parametric sigmoid that is describ ed by Eq.3. The FD is recov ered as the slop e of the curv e in the central p oint at ( x 0 , y 0 ) , according to Eq.4. Based on the box-coun ting approach describ ed in section 3.2, the ‘ds101’ and ‘ds105’ datasets were analyzed for multiple sub jects and v arious choices of smo othing kernel size (see sections 4.2 and 4.2, resp ectively , for details). T able 1 presents the FD estimations for the ‘ds101’ dataset; T able 2 presen ts the FD estimations for the ‘ds105’ dataset. The results include the mean v alues and the corresp onding confidence range at the significance level a = 0 . 05 , as well as the standard deviations. Plain v alues corresp ond to trimmed sets excluding the smallest and largest v alue, while v alues in paren theses corresp ond to the non- trimmed (complete) sets. In b oth tables, each cell corresp onds to FD estimation in 18 instances (9x2 for ‘ds101’ and 6x3 for ‘ds105’). 20 5. EXPERIMENTS AND RESUL TS Fig. 10. The log-log plot log ( P C ( r )) versus l og ( 1 / r ) (blue p oints) and the b est-fit parametric sigmoid (red curv e) that reco vers the fractal dimension (‘FDE’) of the ‘ds105’ dataset (smo othed, sm=8mm 3 ). T able 1. FD estimation in the ‘ds101’ dataset for v arious smoothing k ernel (sm) sizes. Plain v alues corresp ond to trimmed sets excluding the smallest and largest v alue, while v alues in paren theses corresp ond to the non-trimmed (complete) sets. Each cell corresp onds to FD estimation in 18 instances (9x2). Confidence range for mean v alue is at the significance lev el a = 0 . 05 . mean ( µ ) conf.range ( µ ± ) stdev ( σ ) (no sm) 61.07 (60.43) 11.73 (12.57) 23.93 (27.20) sm=4mm 3 31.92 (31.59) 5.13 (5.62) 10.48 (12.17) sm=8mm 3 11.27 (17.14) 2.15 (2.49) 4.39 (5.38) T able 2. FD estimation in the ‘ds105’ dataset for v arious smoothing k ernel (sm) sizes. Plain v alues corresp ond to trimmed sets excluding the smallest and largest v alue, while v alues in paren theses corresp ond to the non-trimmed (complete) sets. Each cell corresp onds to FD estimation in 18 instances (6x3). Confidence range for mean v alue is at the significance lev el a = 0 . 05 . mean ( µ ) conf.range ( µ ± ) stdev ( σ ) (no sm) 15.38 (17.38) 2.86 (5.56) 5.83 (12.04) sm=4mm 3 12.79 (13.75) 2.20 (3.31) 4.48 (7.16) sm=8mm 3 10.67 (11.01) 1.59 (1.89) 3.24 (4.09) 21 6. DISCUSSION 6 Discussion The results presented in section 5.2, as well as the ICA unmixing mo dels that w ere presented in section 5.1, v erify that there is indeed a limited num b er of activ ated brain areas during standard cognitiv e pro cesses. Since these activ ations are present simultaneously , they pro vide a hint of ho w man y tasks are ‘running’ in the h uman brain in p ar al lel as part of its every day functionality . In section 5.1, the results from exp eriments with simulated fMRI data illus- trate the basic unmixing problem for brain sensory data, whic h is relev ant not only to fMRI but other modalities to o, e.g. in EEG. The results sho w that ICA can indeed address the unmixing task with mo derate to go o d p erformance, es- p ecially with regard to the signal sources related to well-defined external stimuli (see comp onen t no.7 in Figure 1 and Figure 2). Due to the nature of ICA and its inheren t constraints, not all signal sources can b e correctly iden tified and, hence, the recov ered comp onents do not match the original ones p erfectly; how ever, if the statistical assertions ab out the signal sources are satisfied adequately , the total reconstruction error can b e minimized effectiv ely . F or the simulated fMRI dataset, the total RMSE for the ICA mixture, reconstructing the original sig- nal with all the recov ered (eight) comp onents, is practically zero (see Figure 3). The most imp ortant results in this case are: (a) the n umber of ICA comp onen ts reco vered matches the num b er of true sources used to construct the original mixture and (b) one of the recov ered comp onents closely matc hes (highly corre- lated) with the w ell-defined external stimuli (square-shap ed time course). This is extremely imp ortant in real fMRI exp erimental proto cols, where sp ecific stimuli t yp es are to b e correlated to sp ecific brain areas for constructing ‘global’ brain atlases . The ICA exp eriments with the real fMRI datasets ‘ds101’ and ‘ds105’, de- scrib ed in sec tion 5.1, illustrate the true p erformance of ICA in constructing factorizations for real brain data. Here, the data volume is muc h larger than in the case of simulated data, since the vo xel grid is now 3-D instead of a sin- gle 2-D ‘slice’, while at the same time the inherent statistics are m uch more complex, as exp ected. F rom Figure 5 and Figure 6 it is clear that ICA w orks as exp ected, providing ‘dense’ (non-sparse) unmixing models with satisfactory p erformance; how ever, it is not alwa ys clear what is the nature of eac h of the reco vered components and ho w they can b e interpreted, especially when sp e- cific signal sources are in question other that a pre-defined external stimuli (e.g. scanner drift, electronic noise, head mov ements, respiration, cardiac pulsation, etc). As describ ed in sections 3.3 and 3.1, in the case of real fMRI datasets the ICA factorization is only approximate (RMSE is nev er zero) and the minimum reconstruction error is achiev ed only when using the maximum allo wable n umber of components - whic h, in turn, is ICA-limited by the num b er of time points a v ailable (i.e., t in matrix Y ∈ R t × n ). In other w ords, a ‘p erfect’ unmixing mo del in real brain data requires the largest p ossible num b er of comp onen ts to b e retrieved. On the other hand, from Figure 7 and Figure 8 it is clear that the reconstruction error drops sharply when the num b er of used comp onen ts is 22 6. DISCUSSION m uch low er than this upp er limit. F or the ‘ds101’ dataset, this num b er seems to b e somewhere in 25 < p < 50 (see Figure 7), while for the ‘ds105’ dataset it is p ' 50 (see Figure 8). In b oth cases, the non-smo othed v ariants of the datasets w ere used, hence there is no loss of fine-detail activ ations and these estimations can b e considered as realistic and consisten t. With regard to the fractal analysis on the sim ulated fMRI data, results in section 5.2 illustrate the robustness and consistency of this metho d. Figure 9 presen ts the log-log plot used to estimate the FD in this case, i.e., the in trinsic dimension of the dataset, which is calculated as 3.83 ( ± 1.45%). This v alue is consisten t with the results of other studies using the same dataset with sparsity- a ware realizations of factorization mo dels [29], where the estimated sparsity is clearly low er (6 or less) than the num b er of signal sources used in the original mixture. F urthermore, Figure 9 sho ws the robustness of the method, with the use of a parametric sigmoid function and T ukey windo w, even when the log-log plot do es not pro vide a clear hin t for the selection of the linear part from where the slop e should be extracted. F or the real fMRI datasets, T ables 1 and 2 present the detailed estimations of the FD for smo othed and non-smo othed v ariants. Specifically , T able 1 sho ws the mean FD v alues for the ‘ds101’ dataset, including the confidence range and the standard deviation. It is clear that, even in the non-smo othed ‘noisy’ v ariant, the in trinsic dimension of the space spanned b y the v o xel data is muc h lo wer ( 48 < D < 63 ) than the dimension of the em b edding space ( 27 K < n < 39 K ). F urthermore, the v alue of FD b ecomes smaller, as exp ected, when smo othing is applied to the data prior to the fractal analysis process. This pro ves that the metho d is consisten t in terms of following the decreasing ‘complexity’ of the dataset, as well as the fact that smo othing the fMRI vo xel data can enhance the qualit y of the most imp ortan t information conten t (ma jor brain activity areas), with a p ossible loss in fine details and/or low-lev el activ ations. Hence, smo othing options in fMRI should b e carefully selected in relation to the sp ecifications of eac h task, i.e., sensitivity versus sp ecificity requirements. Similar comments are v alid for the ‘ds105’, according to the results in T able 2. In the non-smo othed ‘noisy’ v ariant, the intrinsic dimension of the space spanned b y the v oxel data is, again, muc h low er ( 12 < D < 19 ) than the dimension of the em b edding space ( 22 K < n < 48 K ) and it b ecomes even smaller, as exp ected, when smo othing is applied to the data prior to the fractal analysis pro cess. Figure 10 shows the log-log plot used to calculate the FD v alue for the smo othed v ariant (sm=8mm 3 ) of the dataset, where it is clear that the prop osed fractal analysis metho d provides a very reliable estimation. F urthermore, it shows a muc h b etter fit in the sigmoid curve, which means that the b ox-coun ting metho d (see section 3.2) b ecomes more reliable, as expected, w hen the fMRI data are smo othed. As it was men tioned earlier, this study fo cuses on the estimation of the level of p ar al lelism when the h uman brain is p erforming complex cognitive tasks. In some very abstract sense, this is not muc h different than trying to recov er the (minim um) num b er of actual ‘cpu cores’ required to ‘run’ all the active cognitive 23 7. CONCLUSION tasks that are registered in the entire 3-D brain volume while performing a t ypical fMRI exp erimen tal proto col that includes visuo-motor tasks. It is v ery interesting to see that the real fMRI dataset ‘ds101’, which corre- sp onds to a visuo-motor task, pro duces muc h higher estimated FD v alues than the corresp onding FD v alues for the ‘ds105’, which is a muc h simpler visual recognition-only task. This means that, as expected, in the second task there is a muc h low er num b er of distinct activ ated brain areas, hence muc h few er in- dep enden t cognitive tasks inv olved, when no motor resp onse is required by the exp erimen tal proto col. This do es not mean that the total volume of brain ac- tiv ation is smaller but rather than fewer functionalit y comp onents (‘sources’) are present in p ar al lel when visual recognition is concerned, rather than when a prop er motor resp onse is required by the sub ject. This is inherently the casual link to ‘cpu cores’, where several pro cesses are enabled to run simultaneously in a digital computer. ICA reconstruction plots sho w that when the human brain is concerned, this n umber is not defined as a strict threshold but rather in a con tinuous range; when a sp ecific activ ation level is defined, a corresp onding n umber of ‘brain cores’ can b e ev aluated. Ho wev er, in real fMRI data, this range seems to b e non-linear and suc h a num b er can b e retriev ed at the p oint b eyond whic h adding more comp onents has only marginal impact to the mo deled brain signal (see Figures 7 and 8). In short, it seems that normal brain functionality , such as in typical visual or visuo-motor tasks, inv olves only a limited num b er of indep endent pro cesses that run in parallel. Some of them are related to this specific task, while others corresp ond to basic low-lev el functionality , e.g. respiration. Although it is diffi- cult to correctly iden tify and explain all these comp onents in strictly data-driven approac hes (esp ecially in BSS methods like ICA), the inv estigation of the num- b er of ma jor comp onents, in com bination with non-parametric dimensionalit y reco very metho ds such as dataset fractal analysis, can pro vide very useful hints for dev eloping brain-like technologies and algorithms. Curren t research endeav ors lik e the Human Brain Pro ject (HBP) by EU [18] and Brain Research through Adv ancing Inno v ativ e Neurotechnologies (BRAIN) b y USA [47], as w ell as new innov ative VLSI technologies like ‘T rueNorth’ pro ject b y IBM [37,44], require reliable ev aluations of brain activit y not only in the structural but also in the functional level. A typical v oxel size of 3x3x3.5-5 mm 3 corresp onds to roughly 2.5-4 million neurons of sev eral thousands of synapse in terconnections each, or 1 / 40000 to 1 / 25000 of the total brain volume, while the curren tly state-of-the-art ‘T rueNorth’ chip provides 1 million artificial neurons with only 256 synapses each. Hence, the level of true parallelism in human brain is a design asp ect of paramoun t imp ortance in future pro jects. 7 Conclusion This study presents a purely data-driven approach to the estimation of the lev el of p ar al lelism in h uman brain. Using fMRI as the main modality , the h uman brain activit y was inv estigated through ICA for BSS, as well as dataset fractal 24 7. CONCLUSION analysis for the estimation of the intrinsic (true) dimensionality of fMRI data. In some very abstract sense, this is not muc h different than trying to recov er the (minim um) num b er of actual ‘cpu cores’ required to ‘run’ all the active cognitive tasks that are registered in the entire 3-D brain volume while performing a t ypical fMRI exp erimen tal proto col that includes visual-only or visuo-motor tasks. Analysis of the non-smo othed v ariants of the real fMRI datasets (i.e., no information loss) prov ed that even when p erforming complex visuo-motor tasks, the num b er of independent brain processes are in the order of 50 to 60, while it b ecomes muc h low er when visual recognition tasks (no motor response) is concerned. This means that, in theory , an artificial equiv alent of a brain-like cognitiv e structure may not require a massively parallel arc hitecture at the level of single neurons, but rather a prop erly designed set of limited pro cesses that run in parallel on a muc h low er scale. Hence, although current state-of-the-art VLSI technologies still include very limited features and pro cessing p ow er when compared to the real human brain, the assertion of emplo ying actual parallelism lev el of muc h low er order can provide useful hints to future pro jects. A ckno wledgments The author wishes to thank professor Sergios Theo doridis, Yiannis Kopsinis and Anastassios F ytsilis, colleagues at the Dept. of Informatics & T elecommunications, National Kap o distrian Univ. of Athens (NKUA/UoA), for their collab oration in related pro jects and the useful discussions for w orks-in-progress in fMRI analysis and sparse modeling. References 1. V ahid Abolghasemi, Saideh F erdowsi, and Saeid Sanei. Blind separation of image sources via adaptiv e dictionary learning. IEEE T r ansactions on Image Pro c essing , 21(6):2921–2930, June 2012. 2. V ahid Ab olghasemi, Saideh F erdowsi, and Saeid Sanei. F ast and incoherent dic- tionary learning algorithms with application to fMRI. Signal, Image and Vide o Pr oc essing , F ebruary 2013. 3. T rine Julie Abrahamsen and Lars Kai Hansen. Sparse non-linear denoising: Gen- eralization p erformance and pattern reproducibility in functional MRI. Pattern R ec o gnition L etters , 32(15):2080–2085, No vem b er 2011. 4. B. Abrahao and L. Barb osa. Characterizing datasets using fractal metho ds. Dept. of Computer Science, Universidade F e der al de Minas Ger ais, Belo Horizonte, MG Br azil , 2003. 5. M. Behro ozi, M. R. Daliri, and H. Boy aci. Statistical analysis metho ds for the fmri data. Basic and Clinic al Neur oscienc e , 2(4):67–74, 2011. 6. A. J. Bell and T.J. Sejnowski. An information maximisation approach to blind separation and blind deconv olution. Neur al Computation , 7:1129–1159, 1995. 7. V. D. Calhoun, T. Adali, J. Hansen, L. K.and Larsen, and J. J. P ek ar. Ica of functional mri data: An ov erview. In Pr o c. Int. Symp. Indep. Comp. Analysis (ICA 2003) . 8. V. D. Calhoun, T. Adali, G. D. Pearlson, and J. J. Pek ar. Spatial and temp oral indep enden t comp onen t analysis of functional mri data con taining a pair of task- related wa v eforms. Human Brain Mapping , 13(1):43–53, May 2001. 25 7. CONCLUSION 9. C. C. Chen, J. S. Dap onte, and F ox M. D. F ractal feature analysis and classification in medical imaging. IEEE T r ans. Me d. Im. , 8(2):133–142, 1989. 10. Gang Chen, Ziad S. Saad, Jennifer C. Britton, Daniel S. Pine, and Rob ert W. Cox. Linear mixed-effects modeling approach to FMRI group analysis. Neur oImage , 73:176–190, June 2013. 11. M. K. Chung. Gaussian kernel smo othing. L e ctur e notes , pages 1–10, 2012. 12. M. K. Chung. Statistic al and Computational Metho ds in Br ain Image Anal ysis . CR C Press, USA, 2013. 13. N. Correa, T. Adali, Yi-Ou Li, and V.D. Calhoun. Comparison of blind source separation algorithms for fmri using a new matlab to olb ox: Gift. In A c oustics, Sp ee ch, and Signal Pr o c essing, 2005. Pr o c ee dings. (ICASSP ’05). IEEE Interna- tional Confer enc e on , volume 5, pages v/401–v/404 V ol. 5, Marc h 2005. 14. K. P . Cosgrov e, C. M. Mazure, and J. K. Staley . Evolving kno wledge of sex dif- ferences in brain structure, function, and c hemistry . Biol Psychiat , 62(8):847–855, 2007. 15. ARP (DoE-USA). Artificial retina pro ject (arp) – h ttp://artificialretina.energy .gov/. 16. J. L. Encarnacao, H.-O. P eitgen, G. Sak as, and G. (Eds.) Englert. F r actal Ge ometry and Computer Gr aphics . Springer-V erlag, Berlin, 1992. 17. Erik Barry Erhardt, Sriniv as Rac hakonda, Edward J. Bedric k, Elena A. Allen, T ÃŒla y Adali, and Vince D. Calhoun. Comparison of multi-sub ject ica metho ds for analysis of fmri data. Human Br ain Mapping , 32(12):2075–2095, 2011. 18. HBP (EU). The human brain pro ject (hbp) – http://www.h uman brainpro ject.eu/. 19. H. Georgiou, M. Ma vroforakis, N. Dimitrop oulos, D. Cav ouras, and S. Theo dor- idis. Multi-scaled morphological features for the characterization of mammographic masses using statistical classification schemes. Artificial Intel ligence in Me dicine , 41(1):39–59, 2007. 20. S. J. Hanson, T. Matsuk a, and J. V. Haxb y . Com binatorial co des in ven tral tem- p oral lob e for ob ject recognition: Haxby (2001) revisited: is there a ‘face’ area? Neur oimage , 23(1):156–166, 2004. 21. F. J. Harris. On the use of windows for harmonic analysis with the discrete fourier transform. Pro c.IEEE. , 66:66–67, 1978. 22. J. Haxby , M. Gobbini, M. F urey , A. Ishai, J. Schouten, and P . Pietrini. 23. J. V. Haxby , M. I. Gobbini, M. L. F urey , A. Ishai, J. L. Schouten, and P . Pietrini. Distributed and ov erlapping representations of faces and ob jects in ven tral tempo- ral cortex. Scienc e , 293(5539):2425–2430, 2001. 24. A. Hyv ärinen. F ast and robust fixed-p oint algorithms for indep endent comp onent analysis. Neural Networks , 10:626–634, 1999. 25. A. Hyv arinen, J. Karh unen, and E. Oja. Indep endent Component Analysis . Wiley- In terscience, 2001. 26. A. Hyv ärinen and E. Oja. A fast fixed-p oint algorithm for indep enden t comp onen t analysis. Neural Computation , 9(7):1483–1492, 1997. 27. Andreas Jung. An in tro duction to a new data analysis tool: Indep endent comp o- nen t analysis. In Pr o c e e dings of W orkshop GK" Nonline arity"-R e gensbur g , 2001. 28. A. M. Kelly and Milham M. P . ds101: Nyu simon task – op enfmri.org (dataset), 2011. 29. Y annis K opsinis, Harris Georgiou, and Sergios Theodoridis. fmri unmixing via prop erly adjusted dictionary learning. In 22th Eur op e an Signal Pr o c essing Confer- enc e (EUSIPCO 2014) . 30. N. Lazar. The Statistic al Analysis of F unctional MRI Data . Springer, 2008. 26 7. CONCLUSION 31. K. Lee and J. C. Y e. A data-driven fmri analysis using k-svd sparse dictionary learning. International So ciety of Magnetic R esonanc e in me dicine (ISMRM) , 2010. 32. Kang joo Lee, Sungho T ak, and Jong Chul Y e. A data-driven sparse GLM for fMRI analysis using sparse dictionary learning with MDL criterion. IEEE T r ansactions on Me dic al Imaging , 30(5):1076–1089, Ma y 2011. 33. Martin A. Lindquist. The statistical analysis of fMRI data. Statistical Scienc e , 23(4):439–464, Nov em b er 2008. 34. P . R. Massopust. F ractal F unctions, F r actal Surfac es and W avelets . Academic Press, San Diego, 1994. 35. M. Mavroforakis, H. Georgiou, N. Dimitrop oulos, D. Cav ouras, and S. Theo doridis. Significance analysis of qualitative mammographic features using linear classifiers, neural net works and support v ector mac hines. Eur op e an Journal of R adiolo gy , 54:80–89, 2005. 36. M. Mavroforakis, H. Georgiou, N. Dimitrop oulos, D. Cav ouras, and S. Theo doridis. Mammographic masses c haracterization based on localized texture and data set fractal analysis using linear, neural and support vector machine classifiers. Artificial Intel ligenc e in Me dicine , 37(2):145–162, 2006. 37. P aul A. Merolla, John V. Arthur, Ro drigo Alv arez-Icaza, Andrew S. Cassidy , Jun Sa wada, Filipp Akop yan, Bryan L. Jackson, Nabil Imam, Chen Guo, Y utak a Nak a- m ura, Bernard Brezzo, Iv an V o, Stev en K. Esser, Rathinakumar Appusw amy , Brian T aba, Arnon Amir, Myron D. Flickner, William P . Risk, Ra jit Manohar, and Dhar- mendra S. Mo dha. A million spiking-neuron in tegrated circuit with a scalable comm unication netw ork and in terface. Scienc e , 345(6197):668–673, 2014. 38. MLSP-Lab. Mac hine learning for signal pro cessing lab oratory – h ttp://mlsp.umbc.edu/. 39. A. W. Moore and M. S. Lee. Efficient algorithm for minimizing cross v alidation error. In Pr o c e e dings of the 11th.International Conferenc e on Machine L e arning (ML-94), 10-13 July 1994, New Brunswick, New Jersey . 40. National Institute of Mental Health (NIMH) and National Institute of Neurological Disorders Stroke (NINDS). Neuroimaging informatics technology initiative (nifti) – http://nifti.nimh.nih.go v/bac kground. 41. A. J. O’T o ole, F. Jiang, H. Abdi, and J. V. Haxby . Partially distributed rep- resen tations of ob jects and faces in v entral temp oral cortex. J Co gn Neur osci. , 17(4):580–590, 2005. 42. A. Paren t and M. B. Carp enter. Carp enter’s Human Neur o anatomy . Williams- Wilkins, 1995. 43. S. Pierre and J. F. Riv est. On the v alidity of fractal dimension measurements in image analysis. J. Vis. Com. Im. Pr o c. , 7(3):217–229, 1996. 44. A. Sebastian. Ibm cracks op en a new era of computing with brain- lik e chip: 4096 cores, 1 million neurons, 5.4 billion transistors – h ttp://www.extremetech.com/extreme/187612-ibm-crac ks-op en-a-new-era- of-computing-with-brain-lik e-chip-4096-cores-1-million-neurons-5-4-billion- transistors. 45. S. M. Smith. Overview of fmri analysis. The British Journal of R adiolo gy , 77:S167– S175, 2004. 46. C. Jr. T raina, A. J. T raina, L. W u, and C. F aloutsos. F ast feature selection using fractal dimension. In XV Br azilian Datab ase Symp osium 2000, Jo ao Pesso a, P A Br azil . 47. BRAIN (USA). Brain researc h through adv ancing inno v ative neurotechnologies (brain) – http://www.braininitiativ e.nih.gov/. 27