Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Novel region based fMRI analysis using invariant moment descriptors Ng, Bernard 2007

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2007-0532.pdf [ 6.95MB ]
Metadata
JSON: 831-1.0100697.json
JSON-LD: 831-1.0100697-ld.json
RDF/XML (Pretty): 831-1.0100697-rdf.xml
RDF/JSON: 831-1.0100697-rdf.json
Turtle: 831-1.0100697-turtle.txt
N-Triples: 831-1.0100697-rdf-ntriples.txt
Original Record: 831-1.0100697-source.json
Full Text
831-1.0100697-fulltext.txt
Citation
831-1.0100697.ris

Full Text

Novel Region Based fMRI Analysis Using Invariant Moment Descriptors by Bernard Ng B.A.Sc, Simon Eraser University, 2005  A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E D E G R E E O F  M a s t e r o f A p p l i e d Science in. T H E F A C U L T Y O F G R A D U A T E STUDIES (Electrical and Computer Engineering)  The University of British Columbia August 2007 © Bernard Ng, 2007  A b s t r a c t  A new functional magnetic resonance imaging (fMRI) analysis approach based on invariant spatial descriptors was developed for quantitative characterization of brain activation patterns within a region of interest (ROI). In particular, the feasibility of using three dimensional moment invariants (3DMIs) to perform such spatial characterization was examined. The use of spatial descriptors is particularly novel in the field of ROTbased fMRI analysis, since up till now, only magnitude-based features were traditionally employed, which neglect the information encoded by voxel locations within an ROI. The invariance properties of the proposed descriptors to similarity transformations account for inter-subject variability in brain size and subject's orientation within the M R scanner, thus allowing for spatial distributions of activation statistics to be meaningfully compared across subjects. Enhanced sensitivity in detecting task-related activation differences as compared to traditional magnitude-based methods was demonstrated with real fMRI data. To handle the issue of feature selection, a modified linear discriminant analysis (LDA) procedure that incorporates leave-one-out cross-validation was developed. Also, methods to deal with the two main issues in ROI-based fMRI group analysis, namely errors in ROI delineation and inclusions of voxels falsely deemed active, were proposed. One method involves remapping the coordinate space with a Gaussian  ii  function, which in effect de-emphasizes voxels near the R O I boundary, thus also accounts for inter-subject variability in brain shapes. The other method detects outlier voxels that exhibit disproportional influence on the proposed invariant spatial descriptors, and deweights or removes those voxels accordingly. Testing these processing methods on real fMRI data showed further increase in discriminability of task-related activation differences compared to the original 3DMIs alone. To fully exploit the spatio-temporal structure inherent in fMRI data, we extended our spatial characterization approach into the spatio-temporal domain. We showed, for the very first time, that the modulation of the spatial distribution of B O L D signals does in fact correlate with the stimulus, and provides greater sensitivity in detecting activated ROIs and in discriminating task-related differences as compared to traditional mean intensity-based methods. The neuroscience implications of our findings are substantial, and might hence provide brain researchers and clinicians a new promising direction to explore.  iii  Contents Abstract  ii  Contents  iv  List of Tables  ix  List of Figures  x  Acknowledgements  xiv  Statement of Co-authorship  xv  1 Introduction 1.1  Problem Statement and Motivation  2  1.2  Thesis Objective  5  1.3  Magnetic Resonance Imaging  5  1.3.1  Nuclear Magnetic Resonance Phenomenon  6  1.3.2  Magnetic Resonance Imaging Principle  8  1.4  /  1  Functional Magnetic Resonance Imaging  10  1.4.1  Hemodynamic Response  11  1.4.2  Noise Issues in fMRI  13 iv  1.5  1.6  1.7  1.8  fMRI Preprocessing  14  1.5.1  Spatial Realignment and Motion Correction  14  1.5.2  Slice Timing Correction  16  1.5.3  Temporal Autocorrelation Correction  16  1.5.4  Spatial Normalization  17  1.5.5  Spatial Smoothing  19  fMRI Analysis Approaches  20  1.6.1  Exploratory Approach  20  1.6.2  Hypothesis-based Approach  22  Activation Detection using Inferential Statistics  23  1.7.1  Voxel-based Approach  24  1.7.2  ROI-based Approach  26  Contributions to ROI-based fMRI Analysis  Bibliography  27  31  2 fMRI Analysis with 3D Moment Invariants  43  2.1  Introduction  43  2.2  Imaging and Data Pre-processing  45  2.2.1  Study Subjects and Experimental Conditions  45  2.2.2  fMRI Data Acquisition  45  2.2.3  fMRI Pre-Processing  46  2.3  Methods  46  2.3.1  Framework of Proposed Method  46  2.3.2  3D Moment Invariants  47  2.3.3  Linear Discriminant Analysis (LDA)  51  v  2.3.4  Permutation Test  52  2.3.5  Visualization of Activation Patterns  53  2.3.6  Comparison of ROI Analysis Methods  53  2.4 Results and Discussion  55  2.4.1  Supplementary Motor Areas  57  2.4.2  Cerebellum  57  2.4.3  Primary Motor Cotex  58  2.4.4  Prefrontal Cortex  58  2.4.5  Caudate  58  2.5 Conclusion  59  Bibliography 3  60  Feature Selection w i t h M o d i f i e d L D A  64  3.1 Introduction  64  3.2 Materials  69  3.2.1  Study Subjects and Experimental Conditions  69  3.2.2  Data Acquisition  69  3.2.3  fMRI Pre-Processing and Statistical Map Generation  70  3.3 Methods  71  3.3.1  Proposed Framework for Spatial fMRI Analysis  71  3.3.2  ROI Spatial Feature Extraction: 3D Moment Invariants . . .  72  3.3.3  Feature Winsorization  77  3.3.4  Feature Selection: Modified LDA  79  3.3.5  Permutation Test  81  3.3.6  Comparison of ROI Features  82  vi  3.4  Results  83  3.5  Discussion  84  3.6  3.5.1  Neurobiological Relevance of Task-specific ROI Differences  3.5.2  Spatial Analysis of fMRI Data  Conclusions  .  85 86 87  Bibliography  89  4 Joint Spatial Denoising and Active ROI Delineation  95  4.1  Introduction  95  4.2  Data Acquisition  97  4.2.1  Study Subjects and Experimental Conditions  97  4.2.2  Functional M R Imaging  97  4.2.3  Preprocessing and Activation Statistics Generation  98  4.3  Methods  99  4.3.1  ROI Feature Extraction  99  4.3.2  Voxel Influence Map Generation  100  4.3.3  Outlier Voxel Detection  101  4.3.4  Outliner Voxel De-weighting  101  4.3.5  Activation Discriminant Analysis  103  4.3.6  Comparison of Spatial Denoising Methods  103  4.4  Results and Discussion  103  4.5  Conclusion  105 107  Bibliography  vii  5 Temporal Dynamics of Spatial Distributions in fMRI BOLD Signals  109  5.1  Introduction  109  5.2  Data Acquisition and Preprocessing  112  5.2.1  fMRI Data Acquisition  113  5.2.2  fMRI Preprocessing  113  5.3  Methods  114  5.3.1  Feature Time Course Extraction  114  5.3.2  Activation Detection  115  5.3.3  Activation Pattern Discrimination  116  5.4  Results and Discussion  116  5.5  Conclusion  119  Bibliography  120  6 Conclusions and Future Work  122  6.1  Spatial Analysis of Activation Statistics Maps  122  6.2  Spatio-temporal Analysis of BOLD fMRI Signal  124  Bibliography  126  viii  List of Tables 2.1  p-values of Activation Differences During External versus Internal Task 56  3.1  p-values of Activation Differnces During External versus Internal Tasks based on Positive T-statistics  4.1  83  T-values of Activation Map Feature Differences ( A M F D ) Comparing Fast versus Slow Frequency  104  5.1  p-values of ROI Activation  118  5.2  p-values of Activation Differences Comparing Fast versus Slow Frequencies  119  ix  L i s t of Figures  1.1  Regions of interest (ROIs) of motor related brain regions examined in this thesis  1.2  Nuclear Magnetic Resonance,  (a) Orientation of hydrogen atoms  when no magnetic field is present, (b) Orientation of hydrogen atoms when a magnetic field, BQ, is applied, (c) Resulting net magnetic moment, M, in the direction of Bo. (d) Response of M upon excitation by a radio frequency (RF) pulse 1.3  Magnetic Resonance Imaging Principle, (a) Orientation of magnetic moments (black arrows) when a magnetic field, BQ, is applied in the z-direction. (b) Slice selection with an R F pulse. All moments precess at the same frequency, (c) Differences in the rate at which the moments precess when a magnetic gradient is applied in the y-direction. More arrows on the circle indicate faster precession,  (d) Resulting  phase differences between the moments in the y-direction when magnetic gradient is removed, (e) Different precession frequencies when a magnetic gradient is applied in the x-direction.  x  1.4  fMRI Analysis Overview, (a) A series of brain images acquired over time, (b) Intensity changes of a voxel over time. Typically, the intensity time course is contaminated by various noise sources such as low frequency drift, high frequency noise, and motion artefacts, (c) Pre-processed time course compared to an expected response  1.5  Hemodynamic Response resulting from changes in the concentration of de-oxygenated hemoglobin upon stimulus  1.6  11  12  General Linear Model. The matrix Y contains a voxel intensity time course along each column, and the matrix X contains the regressor(s) such as the expected response. Each element in 8 represents the correspondence between a voxel time course and the expected response, and e is the residual matrix  23  2.1  Framework of proposed method  48  2.2  Methods used to compare activation statistics across ROIs. The proposed method captures the rich shape information of the activated regions using 3D moment invariants. Note that irrespective of the orientation of the activation (as demonstrated by the 2 different perspectives shown) the same 3D Moment Invariants will be generated. Conventional methods (right) collapse the spatial information into a single metric such as the mean thresholded T-statistic (upper right) or percentage of activated voxels (lower right), ignoring the texture of the activation, and potentially reducing discriminability across subjects  54  xi  3.1  Simplified synthetic example of ROI activation maps illustrating the underlying idea of the proposed method.  Assuming 100 voxels are  within the ROI and only 3 voxels are active in both conditions. Mean activation statistics = 6/100 and percentage of activation = 3% are the same for both conditions A and B . However, these two tasks are spatially encoded quite differently. This difference in the spatial distribution of activation statistics could be characterized using 3D moment invariants as proposed in this thesis to better discriminate changes in activation patterns between the two conditions  67  3.2  Framework of the proposed method  71  3.3  Spatial feature selection (see description in Section 3.3.4). With the 3DMIs ordered based on their relative discriminability, 5 appears to be the optimal number of features for the dataset examined in this study. When more than five 3DMIs are used, discriminability is reduced due to overfitting  84  4.1  Framework of proposed spatial denoising method  99  4.2  Left thalamus denoised with J\.  Shaded region = outlier voxels  (~20% of the ROI). Most outlier voxels (~80%) lie on the ROI boundary  5.1  105  Experimental task and stimulus timing, (a) Subjects were required to keep the side of the black ring on the gray path (see text), (b) R = rest, Slow, Med, and Fast = stimulus at 0.25, 0.5, and 0.75 Hz, respectively. Each block is 20 s in duration  xii  112  5.2 Feature time course segmentation. The box-car curve corresponds to timing of the stimulus delayed by 4 s. The solid line is a sample feature time course (spatial variance, J\(t), of the left Ml averaged over subjects with its temporal mean removed and divided by its standard deviation). The dotted lines show how the feature time courses are parsed into 6 segments. Slow, Med, and Fast correspond to the task frequencies of 0.25, 0.5, and 0.75 Hz, respectively  xiii  117  Acknowledgements I would like to thank my supervisors Dr. Rafeef Abugharbieh and Dr. Martin J . McKeown for their guidance. I would also like to thank Samantha J . Palmer, Drew Smith, and Cara Slagle for collecting the data used in this thesis. Furthermore, I would like to thank my friends inside and outside BiSICL lab for their support; and last but not least, I would like to thank my parents for their encouragements.  BERNARD N G  The University of British Columbia August 2007  xiv  Statement of Co-authorship All chapters in this thesis were written by Bernard Ng and edited by his supervisors, Dr. Rafeef Abugharbieh and Dr. Martin J . McKeown. The experimental data were collected by Samantha J . Palmer, Drew Smith, and Cara Slagle. A l l data analysis methods were designed and developed by Bernard Ng under the supervision of his supervisors.  xv  Chapter 1  Introduction Traditional neuroimaging mainly focuses on detecting structural abnormalities such as brain atrophy in multiple sclerosis (MS) patients, to detect, track, and evaluate neurological disease. However, examining only the structural properties of the brain is inadequate in explaining many neurological disorders such as Parkinson's Disease (PD), that are associated with the brain's functional  mechanisms.  To facilitate  a broader understanding of human brain functions, considerable research efforts have been invested in advancing functional neuroimaging during the past decade with promising results gradually translated to clinical settings [1]. In this thesis, we present our contributions to one of the most prevalent functional neuroimaging modalities, namely functional magnetic resonance imaging (fMRI). In particular, we introduce a novel region of interest (ROI) based fMRI analysis approach for detecting activated ROIs and discriminating activation pattern changes.  1  1.1  Problem Statement and Motivation  Based on the Global Burden of Disease (GBD) study undertaken by the World Health Organization, World Bank, and Havard School of Public Health, neurological disorders are ranked as the number one disease burden among other diseases examined including the human immunodeficiency virus (HIV), malignant neoplasm, and digestive diseases [2]. Neurological disorders constitutes 6.3% of the global burden of diseases [2], and is predicted to increase as the global population continues to age. Among the neurodegenerative disorders, Parkinson's disease (PD) is found to be the second most common worldwide, affecting 3,765,000 people [3]. In fact, Canada alone has approximately 100,000 people suffering from P D . Thus, in this thesis, an emphasis is placed on characterizing the brain activity of P D patients [4], [5] as compared to healthy subjects [6], [7], [8], [9]. Nonetheless, the generality of our approach permits other neurological diseases to be examined using the same methods presented. P D is a chronic, progressive neurodegenerative disease with primary symptoms such as tremor, rigidity, bradykinesia, and postural instability. These symptoms greatly disrupt the daily activities of P D patients. For instance, even something as simple as eating is affected due to tremoring hands. Also, P D can result in depression, slowed reaction time, difficulty in swallowing, and sleep disturbance among many other complications.  Most motor-related symptoms are believed to  arise from premature death of dopamine-secreting cells in the substantia nigra pars compacta (SNpc) - a brain region involved in motor control. However, the exact cause of premature cell death is unknown, though some studies suggest abnormal accumulation of the protein alpha-synuclein [10]. Currently, the primary clinical tool for diagnosing P D is the Unified Parkin2  son's Disease R a t i n g Scale ( U P D R S ) [11], w h i c h is based o n e x a m i n i n g patients' s y m p t o m s severity.  T h i s measure, however, is insufficient i n fully detecting and  diagnosing P D , p a r t l y due to the brain's compensatory mechanisms, w h i c h obscure the effects of disease.  A l s o , the s y m p t o m s of P D are often confused w i t h the ef-  fects of n o r m a l aging. T o enable better assessment and t r a c k i n g of P D , alternative markers that directly target the source (i.e. the brain) is required. O n e possible biomarker is the functional property of the b r a i n , w h i c h was examined i n this thesis. In p a r t i c u l a r , the b r a i n a c t i v i t y i n motor related regions was studied (see F i g . 1.1).  Supplementary  Prefrontal  F i g u r e 1.1: Regions of interest (ROIs) of motor related b r a i n regions examined i n this thesis.  3  Prior electrophysiological studies have shown that both cortico-cortical [12] and cortico-subcortical coupling [13] are impaired in P D patients, and are modulated by medication [12], [13]. Also, studies on monkeys with P D showed overactivity of the globus pallidus even before symptoms onset [14]. Furthermore, reduced brain activity was found in the supplementary motor area (SMA), dorsolateral prefrontal cortex ( D L P F C ) , and anterior cingulate during voluntary movements for patients with more advanced disease symptoms [15]. A l l these preliminary findings confirm the importance of developing sensitive tools to facilitate characterization of brain functions at various stages of disease. One non-invasive imaging modality that is particularly suitable for studying brain function and permits both spatial and temporal characterization is fMRI. The spatial aspect of fMRI is particularly important, as suggested by McKeown et al. in [16], where brain activation of P D patients was found to focus (i.e. become more localized within a particular brain region) in certain cortical areas upon medication. Thus, extending traditional magnitude-based analysis to the spatial and spatio-temporal domain provides a more complete characterization of brain functions and effects of disease.  Such characterization can be transferred to a wide  range of clinical applications ranging from disease diagnosis to drug therapy evaluation. For example, by examining at which dosage levels does brain activity change significantly, optimal drug dosage can be inferred. Also, earlier P D detection might be feasible by examining how brain function changes with disease progression. This enables clinicians to develop more effective therapies to better treat P D patients at an earlier stage.  4  1.2  Thesis Objective  Traditionally, only magnitude-based features such as mean voxel statistics and percentage of activated voxels [17], are used in ROI-based fMRI analysis.  This lim-  itation on features is mainly due to inter-subject variability, where nonfunctional related differences in brain shapes, brain sizes, and subject's orientation in the scanner will be amplified if spatial information is naively incorporated. However, prior neuroscience studies suggest that different tasks might be spatially encoded by voxel locations within an ROI [16], [18], [19]. Thus, developing methods to quantitatively characterize the spatial distribution of activation statistics can be beneficial. In this thesis, the feasibility of using invariant spatial descriptors to account for inter-subject variability was explored. In particular, three dimensional moment invariants (3DMIs) were examined. Also, processing methods based on the proposed invariant spatial descriptors were developed to account for the two main issues in ROI-based fMRI analysis, namely errors in ROI delineation and inclusions of voxels falsely deemed active. Furthermore, to fully exploit the spatio-temporal structure inherent in fMRI data, we extended our spatial characterization into the temporal domain and explored whether the spatial distribution of the blood oxygenation level-dependent (BOLD) signal itself is modulated in time by task performance.  1.3  Magnetic Resonance Imaging  M R I exploits the nuclear magnetic resonance (NMR) phenomenon inherent to hydrogen atoms to non-invasively and safely image soft tissues in living subjects. Since upon excitation of the tissues, all hydrogen atoms will simultaneously respond, a spatial mapping formed by applying different magnetic gradients at different loca-  5  tions is used to localize each signal source. A brief review of the N M R phenomenon and the M R I principle are provided next.  1.3.1  Nuclear Magnetic Resonance Phenomenon  The human brain is composed of 78% water [20], and each water molecule contains two hydrogen atoms. Each hydrogen atom has a magnetic moment associated with its intrinsic spin. In the absence of an external magnetic field, the orientation of the spin axis of the hydrogen atoms will be random. Thus, no net magnetic moment will be present, as shown in Fig. 1.2(a). When a magnetic field, Bo, is applied, the  (a)  (b)  (c)  (d)  Figure 1.2: Nuclear Magnetic Resonance, (a) Orientation of hydrogen atoms when no magnetic field is present, (b) Orientation of hydrogen atoms when a magnetic field, Bo, is applied, (c) Resulting net magnetic moment, M , in the direction of BQ. (d) Response of M upon excitation by a radio frequency (RF) pulse.  spin axis of the hydrogen atoms will rotate (or precess) around the direction parallel or anti-parallel to Bo, as shown in Fig. 1.2(b). More hydrogen atoms will tend to precess around the parallel direction, since the parallel direction corresponds to a lower energy state, thus resulting.in a net magnetic moment, M (Fig. 1.2(c)). The  6  angular frequency of precession, UQ, is governed by the Larmor relationship:  i/ = 7-fio,  (1-1)  0  where 7 is the gyromagnetic ratio, whose value depends on the nature of the nuclei. For hydrogen, 7 — 42.58 M H z / T . When a radio frequency (RF) pulse matching the precession frequency is applied perpendicular to BQ, M will be flipped towards a direction that is perpendicular to Bo as shown in Fig. 1.2(d). The resulting magnetic moment along the z-direction, M , will be zero while the magnetic moment in the z  xy-plane, M ,  will be equal to M. In this orientation, the hydrogen atoms will be  xy  in an excited state. When the R F pulse is switched off, M  z  will begin to increase  exponentially with a time constant, Tl, also known as spin-lattice relaxation time, as M returns to its equilibrium state (i.e. aligned with Bo)- Tl differs across tissue types. Hence, different tissues will exhibit different M , thereby providing the signal z  contrast required for generating brain images. Another phenomenon that provides signal contrast is the spin-spin interactions of the hydrogen atoms. At equilibrium, hydrogen atoms precess incoherently, but when excited with an R F pulse, the precessions will become in phase, leading to stronger M . xy  When the R F pulse is switched off, the hydrogen atoms will begin  to de-phase exponentially with a time constant, T2, known as spin-spin relaxation time, which also depends on tissue type. field can lead to faster de-phasing, thus M  Also, inhomogeneities of the magnetic xy  in fact decays with a time constant,  T2* that is shorter than T2. This magnetic inhomogeneity-induced contrast is the means by which fMRI images are generated as discussed in Section 1.4.1.  7  1.3.2  Magnetic Resonance Imaging Principle  If only an external magnetic field, BQ, is applied, all hydrogen atoms inside the brain will precess at the same Larmor frequency.  Therefore, exciting the brain  with an R F pulse matched to this Larmor frequency will excite every hydrogen atom. The measured signal thus will be a mixture of signals from multiple regions of the brain. To localize the signal sources, a spatial mapping is created through controlled manipulations of the external magnetic field, as shown in Fig. 1.3. When BQ is applied, the spin axis of all hydrogen atoms in the brain will be aligned in the direction of BQ. If we then apply a magnetic gradient, G , along the z-axis, each z  hydrogen atom will precess at a new frequency governed by:  v = -y(B + G z), 0  (1.2)  z  where z is the z-coordinate of a given group of hydrogen atoms in a small volume of tissue. If we then apply an R F pulse matched to i>, the spin axis of the hydrogen atoms precessing at v will be flipped onto the xy-plane. Thus, a slice of the object will be isolated. After selecting a slice of an object, we still need to isolate each small volume of tissue in the xy-plane. If we apply a magnetic gradient, G , along the y-axis, y  each hydrogen atom will again precess at a new frequency governed by:  v = -y(B + G y), 0  y  (1.3)  where y is the y-coordinate of a given group of hydrogen atoms in the selected slice. Based on (1.3), hydrogen atoms residing at higher gradient will precess faster than those residing at lower gradient. If we then remove the magnetic gradient, all hydrogen atoms will resume to their original precession frequencies with a phase  8  cross sectional view  B e=> 0  A group of neighboring hydrogen atoms In a small volume of tissue  Figure 1.3: Magnetic Resonance Imaging Principle, (a) Orientation of magnetic moments (black arrows) when a magnetic field, Bo, is applied in the z-direction. (b) Slice selection with an R F pulse. A l l moments precess at the same frequency, (c) Differences in the rate at which the moments precess when a magnetic gradient is applied in the y-direction. More arrows on the circle indicate faster precession, (d) Resulting phase differences between the moments in the y-direction when magnetic gradient is removed, (e) Different precession frequencies when a magnetic gradient is applied in the x-direction.  9  difference introduced. Thereby, a mapping between phase to location along the yaxis is generated. If we then apply a gradient along the x-axis, the hydrogen atoms along the higher gradient will again precess at a faster rate. Therefore, by examining the phase and frequency of the signal using the Fourier transform, we can determine the spatial location from which the signal originated. This procedure is repeated at a time interval, TR, for different slices to image a 3D object such as the human brain.  1.4  Functional Magnetic Resonance Imaging  fMRI is based on the assumption that regional changes in blood flow, blood volume or blood oxygenation level during task performance are related to brain activation.  Regional hemodynamic responses can alter the local magnetic fields inside  the brain, which creates image contrast. Among the different contrast techniques, B O L D contrast is the most frequently used [21]. To infer brain activation, multiple TIT-weighted brain volumes are first acquired over time as a subject perform a certain task (Fig. 1.4).  Changes in intensity of a voxel are then tracked over  time to generate an intensity time course. Such time course is then compared to an expected response to infer activation. However, as shown in Fig. 1.4, the raw intensity time course of a voxel can be quite different from that of the expected response. Thus, substantial preprocessing is required to unveil the underlying signal of interest. However, the inherently low signal to noise ratio (SNR) and the complex noise structure of B O L D signal makes brain activation detection rather difficult. In particular, noise sources such as random motion artefacts and location-dependent low frequency drifts are especially hard to handle. The spatio-temporal correlations resulting from these non-functional related confounds can greatly complicate fMRf  10  analysis, and are yet to be resolved.  ^ /' \ . /• v / '  i  /; v  Expected Response  r y A \,  Activated?  0  50  100  150 timers ec)  200  250  Figure 1.4: fMRI Analysis Overview, (a) A series of brain images acquired over time, (b) Intensity changes of a voxel over time. Typically, the intensity time course is contaminated by various noise sources such as low frequency drift, high frequency noise, and motion artefacts, (c) Pre-processed time course compared to an expected response.  1.4.1  Hemodynamic Response  B O L D contrast techniques exploit changes in local magnetic field inhomogenities arising from task-induced modulations of blood oxygenation level, to infer brain activation. The underlying assumption is that the measured B O L D signal is approximately proportional to the amount of local neuronal firing [22]. The temporal profile of B O L D signal (i.e. its hemodynamic response) of a small volume of tissue in response to stimulus is shown in Fig. 1.5. During the onset of stimulus, neurons in a given small volume of tissue will begin to fire, leading to oxygen uptake as 11  Figure 1.5: Hemodynamic Response resulting from changes in the concentration of de-oxygenated hemoglobin upon stimulus.  required for energy metabolism. De-oxygenated blood is paramagnetic thus introduces inhomogeneity to its surrounding magnetic field, resulting in a decrease (an initial dip shown in Fig. 1.5) in the B O L D signal as was discussed in Section 1.3.1. In response to this demand for oxygen, an oversupply of oxygenated blood rushes into the given small volume of tissue. Oxygenated blood is diamagnetic which does not directly affect the local magnetic field; but since the amount of de-oxygenated blood is reduced, the local magnetic field will become relatively more homogeneous, leading to an increase in B O L D signal. After stimulus ceases, the supply of oxygenated blood diminishes as the B O L D signal gradually returns back to baseline after an undershoot.  12  1.4.2  Noise Issues in fMRI  The noise structure of B O L D signals is very complex, which greatly complicates brain activation detection. The signal of interest in the raw intensity time courses are often hidden by many confounds. As shown in the sample raw intensity time course in Fig. 1.4, the signal of interest is contaminated by a mean intensity offset, slow frequency drifts, and high frequency noise in addition to motion artefacts which are harder to detect by inspection.  The magnitude of the mean intensity offset,  however, does not correspond to the level of neural activity.  Instead, this offset  arises from slight T2 differences of each voxel, and is rather difficult to estimate due to the unpredictable trends of the low frequency drifts. Low frequency drifts are suspected to be induced by scanner instability, though the physiological state of the subject may also be a cause.  Among the  different scanner-related noise sources such as magnetic field inhomogeneities and thermal noise, low frequency drifts is probably the most problematic to brain activation detection.  The problem arises from the way brain activation is inferred,  where the higher the correlation between a voxel time course and an expected response signal (see Fig. 1.4), the more likely the voxel is activated. However, low frequency drifts introduce extra autocorrelations to the voxel time courses, which in turn increases correlation, leading to false activation detections.  Methods dealing  with autocorrelation in B O L D signal are discussed in Section 1.5. In addition to scanner noise, fMRI is prone to many physiological noise sources such as motion artefacts, cardiac-linked brain pulsation, respiration, and uncontrolled spontaneous neuronal events. Among these noise sources, motion artefacts is the most challenging to deal with [23].  In general, motion can introduce  artefacts into the B O L D signals by two means. First, inhomogeneities of the mag-  13  netic field induce differences in B O L D signal across positions inside the scanner. Thus, if a subject moves during data acquisition, the B O L D signal of a given voxel will change even if no brain activation is present. Secondly, the strength of the measured B O L D signal depends on spin excitation history [24]. Thus, if a previously excited slice moves into the current acquisition plane, the measured B O L D signal will show reduced amplitude since the previously excited spins would have decayed towards their equilibrium states. This effect, however, is only present for throughplane head motions.  Further complicating fMRI analysis is that subjects tend to  move in response to stimulus [25], thus introducing stimulus-correlated changes in the B O L D signals that are not related to brain activation. A discussion of methods for handling motion-induced artefacts is provided in Section 1.5.  1.5  f M R I Preprocessing  Extracting relevant neurological information from fMRf data is very challenging due to its complex noise structure. Numerous preprocessing procedures are required to reduce the effects of noise in subsequent analyses. Typical procedures include motion correction, slice timing correction, and temporal autocorrelation correction. Also, due to inter-subject variability in brain shape, brain size, and subject's orientation in the scanner, to make group statistical inferences, spatial normalization and spatial smoothing are often performed.  1.5.1  Spatial Realignment and Motion Correction  Most existing motion correction methods involve spatially realigning the motioncorrupted fMRI image volumes to a reference volume, often taken as the first volume or the average of the image volumes. The volumes to be registered to the reference 14  volume are referred to as floating volumes. One approach in realigning the volumes is to define fiducial points and align the fiducial points in the floating volumes to that of the reference volume. However, the low spatial resolution of fMRI images and the lack of distinct fiducial points in the brain makes the fiducial point-based approach inferior to intensity-based approaches [26]. The overall intensity-based realignment process involves three major steps. First, six rigid transformation parameters are estimated to rotate and translate a floating volume so that it becomes aligned with the reference volume. Resampling is then performed to interpolate the intensity values of the floating volume at each voxel location of the reference volume. The optimal resampling method is Fourier interpolation [27], which is equivalent to sine interpolation in the spatial domain. However, sine interpolation involves using every voxel in a volume to interpolate each voxel's new intensity value, thus computationally impractical.  Alternative  interpolation methods include trilinear and B-spline interpolations [28]. Following rigid transformation and resampling, a criterion to determine the spatial correspondence between the floating volume and the reference volume is evaluated. Many criteria have been proposed including image difference [23], image ratio uniformity [29], correlation [30], and mutual information (MI) [31], but no census has been reached as to which criteria is optimal. After criterion evaluation, the described procedures including rigid transformation, resampling, and criterion evaluation are iterated until a desired matching criterion is achieved. Spatial realignment corrects for the differences in orientation of the brain volumes, but does not account for motion-induced intensity changes as discussed in Section 1.4.2. To account for signal changes due to magnetic field inhomogeneities and spin excitation history effects, Friston et al. in [24] proposed using an autore-  15  gressive moving average ( A R M A ) model, which incorporates the position of a given voxel in current and previous scans in correcting the signal intensity of a given voxel. Alternatively, Liao et al. in [32] proposed a method called, "Motion-corrected independent component analysis" ( M C I C A ) , which corrects for motion artefacts by determining a linear combination of spatial transformations applied to certain basis motion-corrupted volumes that maximizes a joint entropy-based criterion. Their simulation results demonstrated that M C I C A is robust to activation level, additive noise, and random motions in the reference volume.  1.5.2  Slice Timing Correction  During analysis, all slices in a brain volume are assumed to be collected simultaneously. However in practice, each slice is collected at a slightly different time. This artefact can lead to phase differences between voxel intensity time courses. For instance, assuming brain activation onsets at time t, if an activated voxel resides in a future slice, the first signal measurement will be made at t + dt. This intensity value will be interpreted as measured at time t, thus the voxel will appear to have begun responding prior to stimulus onset. To account for this artefact, slice timing correction is often performed through Fourier transform, where the magnitude spectrum of the signals is preserved, but the phase is adjusted accordingly. This procedure is equivalent to sine interpolation in time, which is used in standard fMRI analysis packages such as statistical parametric mapping (SPM) [33].  1.5.3  Temporal Autocorrelation Correction  To deal with temporal autocorrelations in the B O L D signals that are not related to brain activation, two approaches are conventionally used - coloring and whiten-  16  ing [34]. Coloring involves high-pass filtering to remove the majority of the low frequency autocorrelations followed by temporally smoothening the BOLD signals with a Gaussian kernel matched to the hemodynamic response function to increase the signal to noise ratio (SNR) [35]. This procedure imposes an autocorrelation structure onto the BOLD signals and presumes that the imposed autocorrelation overwhelms any other noise-induced autocorrelations. The alternative approach is whitening, which involves regressing out the task-related response from the BOLD signals and estimating the autocorrelations in the residual signals for generating a whitening matrix. This matrix is multiplied by the BOLD signal to minimize residual autocorrelations. Existing methods for estimating residual autocorrelation includes AR(1) [36], AR(1) plus white noise [37], AR(p) [38], A R M A [39], and 1// noise model [40]. In addition to the two described approaches, methods based on component analyses have also been recently explored [41]. Such approach involves decomposing the BOLD signals into signals of interest and confounds, as discussed in Section 1.6.1. Also, methods exploiting the whitening or decorrelating property of wavelets have been proposed [42], [43]. 1.5.4  Spatial Normalization  To make group inference in fMRI studies, inter-subject variability in brain shape, brain size, and subject's orientation in the scanner must first be accounted for. One widely used approach is to co-register all subjects' brain images to a common brain template such as the Talairach space [44]. This co-registration is performed by spatially warping or normalizing the structural brain image of each subject to a template, and then applying the same spatial warping to all functional images  17  of the subject. The advantage of warping the structural image instead of directly warping the functional images is that the structural image has significantly higher spatial resolution. The overall spatial normalization procedure involves two steps.  First, an  optimal 12-parameters affine transform is performed to coarsely register the structural image to a template. To ensure that implausible transformations will not be produced, prior information of the parameter values are often incorporated through a Bayesian framework. The additional 6 parameters (zoom and shear along the 3 coordinate directions) compared to rigid transform as discussed in Section 1.5.1, account for brain shape and size differences between a subject's brain and a template. Local refinements can then be performed through non-linear deformations as modeled by a linear combination of smooth spatial basis functions [45]. Alternatively, instead of normalizing the whole brain, approaches for aligning a subject's brain at the ROI level have recently been proposed. Miller et al. [46] used large deformation diffeomorphic metric mappings ( L D D M M ) for ROI registration and demonstrated increased statistical power over conventional whole brain normalization to the Talairach space. A n alternative approach based on continuous medial representations (cm-rep) of the ROIs has been proposed [47]. Several drawbacks to spatial normalization are worth noting. First, the brain shape and size of a subject could be quite different from that of a template, especially for diseased patients. Therefore, improbable deformations of a subject's brain may result. Also, any interpolation involved in the normalization procedure might introduce false information into subsequent analysis. Secondly, estimating a huge number of parameters is typically required during non-linear deformations (ranging from hundreds to thousands), and poor parameter estimates can easily lead to  18  imperfect registration results [48]. Futhermore, normalization artefacts can cause signals from functionally distinct areas, especially small subcortical structures, to be misaligned [49]. Thus, the functional overlap across subjects will be reduced [50]. The most conventional method to deal with this issue is to spatially smooth the data to increase functional overlap. This method, however, may inappropriately pool responses from functionally dissimilar regions, thus reducing the spatial resolution and degrading important spatial information [51]. An alternative method is to adopt an ROI-based approach, where ROIs are manually segmented and statistical properties of regional activation are examined across subjects, as discussed in Section 1.8. This approach has been shown to offer finer localization and increased sensitivity to task-related effects [51]. 1.5.5  Spatial Smoothing  Spatially smoothing the data is intended to serve several purposes. First, the spatial scale of hemodynamic response is approximately 2 to 5 mm, based on high resolution optical imaging experiments [52]. Therefore, by the matched filter theorem, smoothing the data with a kernel of similar size as the hemodynamic response can increase the SNR. Secondly, the Gaussian random field theory used for multiple comparison corrections as to be discussed in Section 1.7.1, requires estimating the spatial smoothness of the data, which is often assumed to be similar to that of the smoothing kernel [52]. Lastly, spatially smoothing the data helps increase the functional overlaps across subjects for making group inferences as discussed in Section 1.5.5. However, this procedure may pool functionally dissimilar regions together, thus rendering subsequent analysis erroneous as pointed out in Section 1.5.4. An alternative approach for increasing SNR without sacrificing as much spa-  19  tial resolution, is to apply spatial wavelet transforms [53], [54]. This approach involves representing clusters of voxels with a few wavelet coefficients and de-weighting or removing those coefficients that fall under a certain threshold (i.e. soft or hard thresholding) chosen based on the level of Type I errors (voxels falsely deemed activated) the experimenter is willing to accept. The inverse wavelet transform is then applied to the surviving coefficients to estimate a denoised version of the data. Also, a multifiltering approach based on combining spatially smoothened data at multiple scales was proposed with enhanced brain activation detection demonstrated [55].  1.6  fMRI Analysis Approaches  In general, fMRI analysis approaches can be classified into two categories, namely exploratory and hypothesis-driven.  The exploratory approach attempts to iden-  tify and group interesting patterns in the data, whereas hypothesis-based approach examines the correspondence between the data and a hypothesized response to experimental stimulus.  1.6.1  Exploratory Approach  The explorative approach seeks to uncover the underling patterns in the data without the need to hypothesize an expected stimulus response, which is often quite difficult. Using this analysis approach, unforeseen signals of interest can be discovered, in addition to task-related signals and confounds. Frequently used methods include clustering [56], [57], [58], [59], [60], [61], [62], [63], principle component analysis ( P C A ) [64], independent component analysis [41], and variants of these methods. Clustering involves identifying and grouping voxels with similar temporal 20  patterns.  Common techniques including K-means [56], fuzzy K-means [57], and  dynamical cluster analysis [59] require specifying the number of clusters and a similarity metric a priori. However, the correct choice of the number of clusters and a similarity metric is not obvious, and the quality of the extracted clusters is difficult to assess. Thus, even though these clustering techniques are computationally efficient, the generated results might be hard to interpret and justify. An alternative approach is component analysis based method, such as P C A and ICA, which involves estimating a matrix, W, to project the voxel intensity time courses onto a space where task-related components and confounds can be separated: S = W-X,  (1.4)  where X is a N x T matrix with each row containing a voxel intensity time course, and S is a JV x T matrix with each row containing the intensity time course of a decomposed component. N and T are the number of voxels within an ROI (e.g. the whole brain) and number of time points, respectively. To determine which voxel might be activated, components that resemble confounds are first removed from S with only task-related components projected back into the voxel intensity time course domain, i.e. X' = W~ S. Non-zero rows of X' thus corresponding to the X  activated voxels. P C A projects the voxel intensity time courses onto a space where the variance of the data is maximized along the projection directions. This criterion lacks neurological motivation, but can nonetheless generate task-related signals and components that resemble scanner and physiological effects as was demonstrated in a study on the effects of nicotine on sustained attention [65]. A more neurologically motivated approach is ICA, which takes advantage of the sparse nature of neural coding [66], i.e. only a portion of the brain activates 21  at a given instant in time. Conventional I C A methods either generate components that are independent in time (mutually independent rows of S) [67] or independent in space (mutually independent columns of W) [41]. However, seeking independent components in time without any spatial constraints can lead to physically improbable forms of W and vice versa [68]. Instead, Stone et al. in [68] proposed using skewed-spatiotemporal I C A to analyze fMRI data, which maximizes statistical independence in both time and space. Task-related time courses generated with skewedspatiotemporal I C A showed superior correspondence with the expected response as compared to P C A , temporal I C A , and spatial I C A [68].  1.6.2  Hypothesis-based Approach  The hypothesis-based approach examines the correspondence between a voxel intensity time course and a hypothesized response to experimental stimulus to infer whether a voxel is activated. The most widely used hypothesis-based method by far is the general linear model ( G L M ) as popularized by the software, S P M [33]: Each column of Y contains a voxel's intensity time course, and each column of the design matrix X contains a regressor, which can be the time courses of the signal of interest or the time course of a confound. Only a boxcar time-locked to stimulus and convolved with the hemodynamic response is illustrated in the example in Fig. 1.6, assuming that motion artefacts and temporal autocorrelations were removed during the preprocessing stage. If confounds were not removed prior to analysis, other regressors can be horizontally concatenated to X to account for those confounds. Each element in 3 represents the correspondence between a voxel's intensity time course (i.e. a column of Y) and a regressor (i.e. a column of X). Each column of  represents the residual of the corresponding voxel. To generate an activation  22  Y  13  f  t  1  <  J  v.  Figure 1.6: General Linear Model. The matrix Y contains a voxel intensity time course along each column, and the matrix X contains the regressor(s) such as the expected response. Each element in (3 represents the correspondence between a voxel time course and the expected response, and e is the residual matrix.  statistics map (or statistical parametric map), each element in (3 is scaled by the variance of the corresponding column of e. A threshold is then set to determine which voxels are activated as to be discussed in Section 1.7.1. We note that hypothesizing X is not trivial.  Addressing this concern, McKeown in [69] proposed a hybrid  approach, where regressors are predicted using spatial I C A .  1.7  A c t i v a t i o n Detection using Inferential Statistics  After an activation statistics map is generated, one can either separately examine whether each voxel is activated (i.e. has statistical value above a certain threshold) or can examine the statistical properties of a group of voxels (ie. an ROI).  23  1.7.1  Voxel-based Approach  T h e voxel-based approach examines each voxel's statistical value separately, thus provides better l o c a l i z a t i o n of b r a i n activation as compared to the R O I - b a s e d approach.  However, to make group inference under such approach, the a c t i v a t i o n  statistics m a p must first be s p a t i a l l y w a r p e d into a c o m m o n space to generate exact correspondence between voxels across subjects. T h i s w a r p i n g procedure c a n introduce artefacts as discussed i n Section 1.5.4.  A l s o , the choice of threshold above  w h i c h a voxel is considered activated is not obvious, though oftentimes a threshold corresponding to an uncorrected p-value of 0.05 is used across a l l voxels, w h i c h is inappropriate due to the m u l t i p l e comparisons problem. W i t h thousands of voxels w i t h i n an activation statistics map, tens or even hundreds of the voxels could easily be falsely deemed active by chance (e.g. 5% of 10,000 voxels = 500 voxels). T o address this problem, most f M R I studies now incorporate some form of m u l t i p l e comparisons correction such as B o n f e r r o n i correction or methods based on r a n d o m field theory [70] or false discovery rate ( F D R ) [71]. Bonferroni correction is based on simple p r o b a b i l i t y rules as follows: L e t QFIV be the family-wise T y p e I error rate, w h i c h is defined as the proba b i l i t y of m a k i n g at least one T y p e I error i n a family of tests (i.e. 1 voxel out of N falsely declared activated) when a l l the null hypothesis are true. If we let a be the p r o b a b i l i t y threshold for each test, then the p r o b a b i l i t y of not m a k i n g a T y p e I error is 1 — a. 1 — (1 — a)  N  F o r N tests, the p r o b a b i l i t y w o u l d be (1 — a) . N  Thus, a^iv =  « N• a since a is s m a l l . T h u s , for 10,000 voxels and CXFW — 0.05 for  example, an a of 0.05/10,000 = 0.000005 is required. However, this threshold is too conservative for most cases, since the voxels are spatially correlated. T o account for this spatial correlation, a less stringent threshold based on G a u s s i a n r a n d o m  24  field  (GRF) theory have been proposed [70]. In brief, determining a threshold based on GRF theory involves first estimating the smoothness of an activation statistics map, which is often assumed to be similar to that of the spatial smoothing kernel. The estimated smoothness value is then used to define the size of a voxel cluster and calculate the expected Euler Characteristic (EC), where E C is the number of voxel clusters that remain after a certain activation statistics threshold is applied. As we increase the threshold, the number of clusters that remain will be either one or zero. Recall that the family-wise Type I error is the probability of making one Type f error in a family of tests. Since the expected EC corresponds to the probability of finding only one cluster above a particular threshold by chance, the expected E C is thus approximately equal to the family-wise error. If we now reverse the argument by first choosing a certain C*FW>  w  e  c  a  n  determine the corresponding activation statistics threshold at which  the family-wise error will be less than apwThe described GRF based procedure requires spatially smoothing the data, which might inappropriately pool functionally dissimilar regions together, as mention in Section 1.5.5. An alternative method that is also less stringent than Bonferroni correction uses FDR to determine which voxels are likely to be activated. F D R is defined as the proportion of false positives (voxels falsely deemed active) among those tests for which the null hypothesis is rejected. Bonferroni correction, on the other hand, controls the amount of false positives among all tests whether or not the null hypothesis is rejected. This FDR-based procedure was introduced in [72] and works as follows:  Let q be the maximum F D R between 0 and 1 that a researcher is willing to ac-  25  cept on average. Order the p-values of the N voxels from smallest to largest:  p(l)  ^ ... ^ p(N). For i from 1 to N, voxels with p(i) ^ i • q/ V are declared activated.  The advantage of the F D R approach is that the threshold is adaptive to each subject's data, thus relieves the need to find an arbitrary threshold that works for all subjects. However, spatial correlation is not accounted for, which implies the calculated threshold will still be too conservative depending on the level of spatial correlation present.  1.7.2  ROI-based Approach  Voxel level inferences have the advantage of providing precise localization of activated voxels within a subject. However, to make group inferences under such scheme, both spatial normalization and spatial smoothing are required which in fact reduces spatial resolution [51]. The alternative approach is to drawn inference at the ROI level, which involves segmenting out ROIs for each subject individually, and examining the statistical properties of regional activations across subjects. This approach has been shown to offer finer localization and increased sensitivity to task-related effects [51]. However, any chosen regional statistics must account for inter-subject variability in brain shape, brain size, and subject's orientation in scanner for inference to be meaningful. Towards this end, past studies have used magnitude-based features, such as mean voxel activation statistics and percentage of activated voxels [17], which neglect any spatial information to avoid complications arising from the differences in subject's orientation in the scanner. The spatial information encoded by the location of voxels within an ROI, however, might be an important attribute of brain activations as suggested by several studies. For example, Thickbroom et  26  al. [18] showed that the spatial extent of activation, as opposed to response magnitude, was modulated by different levels of force during a sustained finger flexion task. Also, McKeown et al. [16] showed by applying I C A to electroencephalogram (EEG) data that L-Dopa medication appears to have a 'focussing' effect on the cortical activity of P D patients. Furthermore, Haxby et al. demonstrated that brain response can be inferred from the pattern of activity even at the macroscopic scale [19], analogous to the idea of population coding where brain response is encoded by the overall pattern of activity in a set of (spatially-disparate) neurons [73]. Thus, developing tools to quantitatively characterize the spatial distribution of activation statistics will definitely be fruitful. Also, a major issue in ROI-based fMRI analysis is errors in ROI delineation. Conventionally, ROIs are manually segmented with the aid of anatomical atlases. However, any errors introduced by the segmentation process will increase inter-subject variability. In this thesis, we propose new methods to incorporate spatial information into ROI-based fMRI analysis that simultaneously accounts for ROI segmentation errors, as to be discussed in the next section.  1.8  Contributions to ROI-based fMRI Analysis  The contributions of this thesis are three-folds: 1. We demonstrated for the very first time that spatial information encoded by the location of voxels within an ROI can be meaningfully incorporated into ROI-based fMRI analysis. 2. Processing methods to deal with ROI delineation errors, brain shape intersubject variability, and inclusion of voxels falsely deemed activated were developed.  27  3. We extended the spatial analysis into the temporal domain, and showed that the spatial distribution of B O L D signal itself is in fact modulated in time by task performance. The use of spatial descriptors is particularly novel in the field of ROI-based fMRI analysis, since up till now, only magnitude-based features are considered. The invariant nature of the proposed spatial descriptors to similarity transformations (i.e. translations, rotations, and scaling) accounts for inter-subject variability in brain size and subject's orientation in the scanner, thus enables spatial characterization without resorting to spatial normalization. In [6], we demonstrated that incorporating spatial information using 3DMIs (as described in detail in Chapter 2 and 3) enhances sensitivity in discriminating task-related activation changes.  Our proposed analysis framework also allows one  to systematically compare two sets of brain activations. across subject groups (e.g.  The comparison can be  diseased group versus control group) or within group  where the same set of subjects perform two different tasks or perform the same task before and after medication. Hence, our method can prospectively be translated to a wide range of clinical applications such as disease diagnosis, prognosis, and drug evaluation. In fact, the method proposed in [6] and its extension [7] were used to examine the effects of L-Dopa medication on P D patients [4], [5]. With infinitely many 3DMIs, deciding which ones to use is not trivial. In [7], we developed a modified L D A procedure that incorporates leave-one-out crossvalidation to deal with this feature selection problem. Furthermore, we modified the 3DMIs by remapping the coordinate space with a Gaussian function.  This  remapping stresses voxels near the activation centroid over voxels near the ROI boundary, which in effect accounts for both ROI delineation errors and inter-subject  28  variability in brain shapes. However, all voxels near the ROI boundary are affected regardless of whether they belong to the ROI. To more directly address the problem associated with ROI delineation errors, a processing method based on 3DMIs was proposed [8] that simultaneously refines ROI delineation and spatially denoises the fMRI activation statistics. The underlying idea is that no voxel in isolation should exhibit exceedingly disproportional influence on the 3DMIs. Thus, voxels deemed overly influential are either deweighted or removed. Applying this procedure to real fMRI data demonstrated further increase in sensitivity for discriminating taskrelated activation differences over using 3DMIs alone. The results in [6], [7], and [8] all suggest that incorporating spatial information into ROI-based fMRI analysis is beneficial. These analyses were performed on T-maps where the spatial information is collapsed over time, thus only the timeaveraged spatial patterns of brain activity were examined. To fully exploit the spatio-temporal structure of fMRI data, we also extended our spatial analysis to the spatio-temporal domain to explore whether the spatial distribution of the BOLD signal itself is modulated in time by task performance [9]. The spatial feature time courses generated by calculating the 3DMIs at each time point, were used to infer ROI activation, in addition to discriminating activation differences as in [6], [7], and [8]. Our results demonstrated that the spatial distribution of BOLD signal is in fact modulated by task performance. The neuroscience implication of this finding is substantial, as this is the first time task-related spatial modulation of brain activation is demonstrated. In all past studies, the traditional belief is that only the magnitude of BOLD signal is modulated by task. This finding can potentially be extended to functional connectivity analysis to more meaningfully summarize the temporal response of an ROI as opposed to averaging the intensity time courses of  29  all voxels within an ROI. All our results suggest the importance of incorporating spatial information in ROI-based fMRI analysis. Exploiting the information encoded by voxel locations might hence provide brain researchers a whole new area of research to explore.  30  B i b l i o g r a p h y  [1] P. M . Matthews, G. D. Honey, and E. T. Bullmore.  Applications of fmri  in translational medicine and clinical practice. Nature reviews Neuroscience, 7(9):732-744, Sep 2006. [2] World Health Organization. Neurological Disorders Public Health Challenges.  Geneva:WHO, 2006. [3] World Health Organization. The world health report 1997. Conquering suffering, enriching humanity. Report of the Director General. Geneva:WHO, 1997.  [4] M . J. McKeown, B. Ng, M . M . Lewis, R. Abugharbieh, and X. Huang. Task and hand dominance-specific "focussing" effect of 1-dopa in parkinson's disease (pd) and normal subjects. 10 International Congress of Parkinson's Disease th  and Movement Disorders, Kyoto, Japan, Oct. 28-Nov. 2, 2006 (abstract). [5] S. J. Palmer, M . J. McKeown, B. Ng, and R. Abugharbieh. 3d-texture of bold fmri activation in parkinson's disease reveals recruitment of ipsilateral cerebello-thalamo-cortical loops modulated by 1-dopa and increasing movement frequency. Annual Meeting: Society for Neuroscience 2007, San Diego, California, Nov. 3-7, 2007 (abstract).  31  [6] B. Ng, R. Abugharbieh, X . Huang, and M . J. McKeown. Characterizing fmri activations within regions of interest (rois) using 3d moment invariants. IEEE Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, NY, June 17-18, 2006, pp. 63 (8 pages). [7] B. Ng, R. Abugharbieh, X . Huang, and M . J. McKeown. Characterization of fmri activation maps using invariant 3d moment descriptors. IEEE on Medical Imaging,  Transactions  invited paper, (under revision).  [8] B. Ng, R. Abugarbieh, S. J. Palmer, and M . J. McKeown. Joint spatial denoising and active region of interest delineation in functional magnetic resonance imaging. 29  th  Annual International Confernce of the IEEE Engineering  in Medicine and Biology Society, Lyon, France, Aug. 23-26, 2007 (accepted, 4 pages). [9] B. Ng, R. Abugharbieh, S. J. Palmer, and M . J. McKeown. Characterizing task-related temporal dynamics of spatial activation distributions in fmri bold signals. 10 International Conference on Medical Image Computing and Comth  puter Assisted Intervention, Brisbane, Australia, Oct. 29-Nov. 2, 2007 (accepted, 8 pages). [10] A. A. Cooper, A. D. Gitler, A. Cashikar, C. M . Haynes, K . J. Hill, B. Bhullar, K. Liu, K. Xu, K . E. Strathearn, F. Liu, S. Cao, K. A. Caldwell, G. A. Caldwell, G. Marsischky, R. D. Kolodner, J. Labaer, J. C. Rochet, N . M . Bonini, and S. Lindquist. Alpha-synuclein blocks er-golgi traffic and rabi rescues neuron loss in parkinson's models. Science, 313(5785):324-328, Jul 21 2006. [11] D. J. Gelb, E. Oliver, and S. Gilman. Diagnostic criteria for parkinson disease. Achives of Neurology,  56(1), 1999.  32  [12] P. Silberstein, A. Pogosyan, A. A. Kuhn, G. Hotton, S. Tisch, A. Kupsch, P. Dowsey-Limousin, M. I. Hariz, and P. Brown. Cortico-cortical coupling in parkinson's disease and its modulation by therapy. Brain : a journal of neurology,  128(Pt 6):1277-1291, Jun 2005.  [13] A. L. Whone, R. L. Watts, A. J. Stoessl, M. Davis, S. Reske, C. Nahmias, A. E. Lang, O. Rascol, M. J. Ribeiro, and P. Remy. Slower progression of parkinson's disease with ropinirole versus levodopa: The real-pet study. Annals of Neurology,  54(1):93-101, 2003.  [14] E. Bezard, A. R. Crossman, C. E. Gross, and J. M. Brotchie. Structures outside the basal ganglia may compensate for dopamine loss in the presymptomatic stages of parkinson's disease. The FASEB journal : official publication of the Federation of American Societies for Experimental Biology,  15(6):1092-1094,  Apr 2001. [15] M. Jahanshahi, I. H. Jenkins, R. G. Brown, C. D. Marsden, R. E. Passingham, and D. J. Brooks. Self-initiated versus externally triggered movements, i. an investigation using measurement of regional cerebral blood flow with pet and movement-related potentials in normal and parkinson's disease subjects. Brain : a journal of neurology, 118(4):913-933, Aug 1995.  [16] M. J. McKeown, S. J. Palmer, W. L. Au, R. G. McCaig, R. Saab, and R. AbuGharbieh. Cortical muscle coupling in parkinson's disease (pd) bradykinesia. Journal of neural transmission. Supplementum,  (70)(70):31-40, 2006.  [17] R. T. Constable, P. Skudlarski, E. Mencl, K. R. Pugh, R. K. Fulbright, C. Lacadie, S. E. Shaywitz, and B. A. Shaywitz. Quantifying and comparing region-  33  of-interest activation patterns in functional brain rnr imaging:  methodology  considerations. Magnetic resonance imaging, 16(3):289-300, Apr 1998. [18] G . W . Thickbroom, B. A . Phillips, I. Morris, M . L . Byrnes, and F . L . Mastaglia. Isometric force-related activity in sensorimotor cortex measured with functional mri. Experimental brain research, 121(l):59-64, Jul 1998. [19] J . V . Haxby, M . I. Gobbini, M . L . Furey, A . Ishai, J . L . Schouten, and P. Pietrini. Distributed and overlapping representations of faces and objects in ventral temporal cortex. Science, 293(5539):2425-2430, 2001. [20] H . Mcllwain and H . S. Bachelard. Biochemistry tem. Churchill Livingstone,  and the Central Nervous Sys-  1971.  [21] S. Ogawa, T . M . Lee, A . R. Kay, and D. W . Tank. Brain magnetic resonance imaging with contrast dependent on blood oxygenation.  Proceedings of the  National Academy of Sciences, 87(24):9868-9872, 1990. [22] D. J . Heeger and D. Ress. What does fmri tell us about neuronal activity? Nature Reviews Neuroscience, 3(2): 142-151, 2002. [23] K . J . Friston, P. Jezzard, and R. Turner. Analysis of functional mri time-series. Human brain mapping, 1(2):153—171, 1994. [24] K . J . Friston, S. Williams, R. Howard, R. S. Frackowiak, and R. Turner. Movement-related effects in fmri time-series. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, 35(3):346-355, Mar 1996. [25] J . V . Hajnal, R. Myers, A . Oatridge, J . E . Schwieso, I. R. Young, and G . M . Bydder. Artifacts due to stimulus correlated motion in functional imaging of 34  the brain. Magnetic  Magnetic Resonance  resonance  in medicine  : official journal  in Medicine / Society of Magnetic  of the Society of  Resonance  in  Medicine,  31(3):283-291, Mar 1994. [26] S. C. Strother, J. R. Anderson, X . L. Xu, J. S. Liow, D. C. Bonar, and D. A. Rottenberg. Quantitative comparisons of image registration techniques based on high-resolution mri of the brain.  Journal  of computer  assisted  tomography,  18(6):954-962, Nov-Dec 1994. [27] W. F. Eddy, M . Fitzgerald, and D. C. Noll. Improved image registration by using fourier interpolation.  Magnetic  nal of the Society of Magnetic Resonance  in Medicine,  resonance  Resonance  in medicine  in Medicine  /  : official  Society of  jour-  Magnetic  36(6):923-931, Dec 1996.  [28] T. M . Lehmann, C. Gonner, and K . Spitzer. Survey: interpolation methods in medical image processing.  IEEE  Transactions  on Medical  Imaging,  18(11):1049-1075, 1999. [29] R. P. Woods, S. T. Grafton, C. J. Holmes, S. R. Cherry, and J. C. Mazziotta. Automated image registration: I. general methods and intrasubject, intramodality validation.  Journal  of computer  assisted tomography,  22(1):139-  152, 1998. [30] D. L. Collins, A. C. Evans, C. Holmes, and T. M . Peters. Automatic 3d segmentation of neuro-anatomical structures from mri. Medical Imaging,  pages 139-152, 1995.  35  Information  Processing  in  [31] F. Maes, D. Vandermeulen, and P. Suetens. Comparative evaluation of multiresolution optimization strategies for multimodality image registration by maximization of mutual information. Medical image analysis, 3(4):373-386, 1999. [32] R. Liao, J. L. Krolik, and M . J. McKeown. An information-theoretic criterion for intrasubject alignment of fmri time series: motion corrected independent component analysis. IEEE Transactions  on Medical Imaging,  24(l):29-44, Jan  2005. [33] K . J . Friston. Statistical parametric mapping and other analysis of functional imaging data. In Brain Mapping:  The Methods,  pages 363-385. Academic Press,  1996. [34] K . J. Friston, 0. Josephs, E. Zarahn, A. P. Holmes, S. Rouquette, and J. Poline. To smooth or not to smooth? bias and efficiency in fmri time-series analysis. Neurolmage,  12(2): 196-208, Aug 2000. .  [35] K . J. Friston, A. P. Holmes, J. B. Poline, P. J. Grasby, S. C. Williams, R. S. Frackowiak, and R. Turner. Analysis of fmri time-series revisited.  Neurolmage,  2(l):45-53, Mar 1995. [36] E. Bullmore, M . Brammer, S. C. Williams, S. Rabe-Hesketh, N . Janot, A. David, J. Mellers, R. Howard, and P. Sham. Statistical methods of estimation and inference for functional mr image analysis. Magnetic medicine  : official journal  Society of Magnetic  of the Society of Magnetic  Resonance  in Medicine,  36  Resonance  resonance  in  in Medicine  /  35(2):261-277, Feb 1996.  [37] P. L. Purdon and R. M. Weisskoff. Effect of temporal autocorrelation due to physiological noise and stimulus paradigm on voxel-level false-positive rates in fmri.  Human  brain mapping,  6(4):239-249, 1998.  [38] K. J. Worsley, C. H. Liao, J. Aston, V. Petre, G. H. Duncan, F. Morales, and A. C. Evans. A general statistical analysis for fmri data.  Neurolmage,  15(1):1-15, 2002. [39] J. J. Locascio, P. J. Jennings, C. I. Moore, and S. Corkin. Time series analysis in the time domain and resampling methods for studies of functional magnetic resonance brain imaging.  Human  brain mapping,  5(3): 168-193, 1997.  [40] E. Zarahn, G. K. Aguirre, and M. DEsposito. Empirical analyses of bold fmri statistics.  Neurolmage,  5(3):179-197, 1997.  [41] M. J. Mckeown, S. Makeig, G. G. Brown, T. P. Jung, S. S. Kindermann, A. J. Bell, and T. J. Sejnowski. Analysis of fmri data by blind separation into independent spatial components.  Human  brain mapping,  6(3): 160-188, 1998.  [42] E. Bullmore, C. Long, J. Suckling, J. Fadili, G. Calvert, F. Zelaya, T. A. Carpenter, and M. Brammer. Colored noise and computational inference in neurophysiological(fmri) time series analysis: Resampling methods in time and wavelet domains.  Human  brain mapping,  12(2):61-78, 2001.  [43] E. Bullmore, J. Fadili, M. Breakspear, R. Salvador, J. Suckling, and M. Brammer. Wavelets and statistical analysis of functional magnetic resonance images of the human brain.  Statistical  methods in medical research,  2003.  37  12(5):375-399,  [44] J. Talairach and P. Tournoux. Co-planar stereotaxic atlas of the human brain: an approach to medical cerebral imaging. Stuttgart, Germany:  Thieme,  1988.  [45] K . J. Friston, J. Ashburner, C. D. Frith, J. B. Poline, J. D. Heather, and R. S. J. Frackowiak. Spatial registration and normalization of images. brain mapping,  Human  2(l):l-25, 1995.  [46] M . I. Miller, M . F. Beg, C. Ceritoglu, and C. Stark. Increasing the power of functional maps of the medial temporal lobe by using large deformation diffeomorphic metric mapping. Proceedings of the National Academy of Sciences, 102(27):9685-9690, 2005. [47] P. A. Yushkevich, J. A. Detre, D. Mechanic-Hamilton, M . A. FernndezSeara, K . Z. Tang, A. Hoang, M . Korczykowski, H. Zhang, and J. C. Gee. Hippocampus-specific fmri group activation analysis using the continuous medial representation. Neurolmage,  35(4):1516-1530, 2007.  [48] O. T. Carmichael, H . A. Aizenstein, S. W. Davis, J. T. Becker, P. M . Thompson, C. C. Meltzer, and Y . Liu. Atlas-based hippocampus segmentation in alzheimer's disease and mild cognitive impairment. Neurolmage,  27(4):979-  990, Oct 1 2005. [49] M . Ozcan, U. Baumgartner, G. Vucurevic, P. Stoeter, and R. D. Treede. Spatial resolution of fmri in the human parasylvian cortex: comparison of somatosensory and auditory activation. Neurolmage,  25(3):877-887, Apr 15 2005.  [50] C. E. Stark and Y. Okado. Making memories without trying: medial temporal lobe activity associated with incidental memory formation during recognition.  38  The Journal of neuroscience : the official journal  of the Society for Neuro-  science, 23(17):6748-6753, Jul 30 2003. [51] A . Nieto-Castanon, S. S. Ghosh, J . A . Tourville, and F . H . Guenther. Region of interest based analysis of functional imaging data. Neurolmage, 19(4): 13031316, 2003. [52] R. S. J . Frackowiak, K . J . Friston, C . Frith, R. Dolan, C . J . Price, S. Zeki, Human Brain Function.  J. Ashburner, , and W . D. Penny.  Academic Press,  2003. [53] D. Van De Ville, T . Blu, and M . Unser.  Integrated wavelet processing and  spatial statistical testing of fmri data. Neurolmage, 23(4):1472-1485, Dec 2004. [54] M . J. Fadili and E . T . Bullmore.  A comparative  evaluation of wavelet-  based methods for hypothesis testing of brain activation maps.  Neurolmage,  23(3):1112-1128, 2004. [55] J . B. Poline and B. M . Mazoyer. Enhanced detection in brain activation maps using a multifiltering approach. Journal of cerebral blood flow and metabolism :  official journal  of the International  Society of Cerebral Blood Flow and  Metabolism, 14(4):639-642, Jul 1994. [56] D. Balslev, F . A . Nielsen, S. A . Frutiger, J . J . Sidtis, T . B. Christiansen, C. Svarer, S. C . Strother, D. A . Rottenberg, L . K . Hansen, O. B. Paulson, and I. Law. Cluster analysis of activity-time series in motor learning. brain mapping, 15(3):135-145, Mar 2002.  39  Human  [57] R. Baurngartner, C. Windischberger, and E. Moser. Quantification in functional magnetic resonance imaging: fuzzy clustering vs. correlation analysis. Magnetic resonance imaging, 16(2):115—125, 1998.  [58] P. Filzmoser, R. Baurngartner, and E. Moser. A hierarchical clustering method for analyzing functional mr images. Magnetic resonance imaging, 17(6) :817826, Jul 1999. [59] A. Baune, F. T. Sommer, M. Erb, D. Wildgruber, B. Kardatzki, G. Palm, and W. Grodd. Dynamical cluster analysis of cortical fmri activation. Neurolmage, 9(5):477-489, 1999. [60] D. Cordes, V. Haughton, J. D. Carew, K . Arfanakis, and K. Maravilla. Hierarchical clustering to measure connectivity in fmri resting-state data. Magnetic resonance imaging, 20(4):305-317, May 2002. [61] X . Golay, S. Kollias, G. Stoll, D. Meier, A. Valavanis, and P. Boesiger. A new correlation-based fuzzy logic clustering algorithm for fmri. Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, 40(2):249-260, Aug 1998.  [62] E. Moser, M . Diemling, and R. Baurngartner. Fuzzy clustering of gradient-echo functional mri in the human visual cortex, part ii: quantification. Journal of magnetic resonance imaging : JMRI, 7(6):1102-1108, Nov-Dec 1997.  [63] S. H. Yee and J. H. Gao. Improved detection of time windows of brain responses in fmri using modified temporal clustering analysis. Magnetic resonance imaging, 20(l):17-26, Jan 2002.  40  [64] A. H. Andersen, D. M . Gash, and M . J. Avison. Principal component analysis of the dynamic response measured by fmri: a generalized linear systems framework. Magnetic resonance imaging, 17(6):795—815, Jul 1999.  [65] N . S. Lawrence, T. J. Ross, R. Hoffmann, H. Garavan, and E. A. Stein. Multiple neuronal networks mediate sustained attention. Journal of cognitive neuroscience, 15(7):1028-1038, 2003. [66] H. B. Barlow. Unsupervised learning. Neural Computation, 1:295-311, 1989. [67] S. Makeig, T. P. Jung, A. J. Bell, D. Ghahremani, and T. J. Sejnowski. Blind separation of auditory event-related brain responses into independent components. Proceedings of the National Academy of Sciences, 94(20): 10979-10984,  1997. [68] J. V. Stone, J. Porrill, N. R. Porter, and I. D. Wilkinson. Spatiotemporal independent component analysis of event-related fmri data using skewed probability density functions. Neurolmage, 15(2):407-421, 2002. [69] M. J. McKeown. Detection of consistently task-related activations in fmri data with hybrid independent component analysis. Neurolmage, ll(l):24-35, 2000. [70] K . J. Worsley, A. C. Evans, S. Marrett, and P. Neelin. A three-dimensional statistical analysis for cbf activation studies in human brain. Journal of cerebral blood flow and metabolism : official journal of the International  Society of  Cerebral Blood Flow and Metabolism, 12(6):900-918, 1992.  [71] C.R. Genovese, N.A. Lazar, and T. Nichols. Thresholding of statistical maps in functional neuroimaging using the false discovery rate. Neuroimage, 15(4):870878, 2002. 41  [72] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: A practical and powerful approach to multiple testing. Journal of the Royal Statistical Society.Series B (Methodological),  57(l):289-300, 1995.  [73] P. Dayan and L. F. Abbott. Theoretical neuroscience: computational and mathematical modeling of neural systems.  MIT Press, 2001.  42  Chapter 2  fMRI Analysis with 3D Moment Invariants 2.1  Introduction  Most conventional functional magnetic resonance imaging (fMRI) analysis methods generate statistics to determine the probability that a given voxel is being activated during the performance of the underlying behavioral task(s). To combine fMRf statistics across subjects, there are generally two broad approaches. One involves warping the brains of each individual subject to a common exemplar shape [1]. The other involves drawing the regions of interest (ROIs) of each individual subject, and then examining the statistical properties of activations in the ROIs across subjects. Determining how to summarize the different activation statistics within an ROI is not obvious. Mean value, percentage of activated voxels, and variations of these ' A version of this chapter has been published. B. Ng, R. Abugharbieh, X. Huang, and M . J. McKeown, "Characterizing fMRI Activations within Regions of Interest (ROIs) using 3D Moment Invariants," I E E E Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, N Y , June 17-18, 2006, pp. 63 (8 pages).  43  two metrics are typically used to represent the level or extent of activation inside a ROI. For instance, Allen and Courchesne [2] used the percentage of activated voxels and mean percent signal changes to study the activation differences between subjects with and without autism. Aizenstin et al. [3] used the mean of the intensity time courses of significantly activated voxels and percentage of voxels with negative T-values to compare the activation distributions of young and elderly subjects. Similarly, Bogorodzki et al. [4] used features extracted from mean time-intensity curves to discriminate between different groups of marijuana smokers.  A l l these  discrimination techniques are essentially based on comparing activation statistics amplitudes, and will be referred to as intensity-based methods for the remainder of this chapter. While this allows the results from different subjects to be combined in an uncomplicated manner, the spatial information related to the location of the activated voxels is ignored. This spatial information, however, could provide more discriminative power over traditional approaches. Thus, we propose a method for characterizing the distribution of fMRI activation statistics of each voxel inside an ROI using three-dimensional (3D) shape descriptors based on 3D moment invariants. The invariant nature of any shape descriptor is critical, as it would not be desirable to have any shape information that was dependent upon the particular coordinate system used.  A recent study by Mangin et al.  [5] demonstrated the  discriminative power of 3D moment invariants on structural data. In this thesis, we proposed extending the application of these invariants to functional data. In this chapter, we show that applying the proposed method to fMRI data recorded from healthy subjects performing externally and internally guided tasks reveals significant handedness-related and task-related activation differences in the supplementary motor areas, cerebellum, primary motor cortex, prefrontal cortex, and caudate. In  44.  contrast, many of these activation differences were not detected using conventional intensity-based methods.  2.2  Imaging and Data Pre-processing  2.2.1  Study Subjects and Experimental Conditions  In this study, fMRI data was collected from 8 healthy subjects. S M A , cerebellum, primary motor cortex, prefrontal cortex, and caudate were chosen as the regions of interest.  Each patient was asked to perform two motor tasks in 20 s blocks first  with their right hand and then with their left hand:  A. In the first task, patients were externally guided, where they had to follow a finger tapping sequence shown on a screen. B . In the second task, patients were internally guided, where they had to perform a finger tapping sequence that they had seen on a screen earlier.  2.2.2  fMRI Data Acquisition  The fMRI data were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangern, Germany) with a birdcage type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The participant's head was positioned along the canthomeatal line. Foam padding was used to limit head motion within the coil. High-resolution T l weighted anatomical images (3D S P G R , TR=14ms, TE=7700 ms, flip a n g l e = 2 5 ° , voxel dimensions 1.0 m m x l . O m m x l . O mm, 176x256 voxels, 160 slices were acquired for co-registration and normalization of functional images.  A total of 49 co-planar functional images were acquired using a gradient  45  echoplanar sequence (TR=3000 ms, T E = 3 0 ms, flip a n g l e = 8 0 ° , N E X = 1 , voxel dimensions 3.0 mmx3.0 mmx3.0 mm, imaging matrix 64x64 voxels). Each functional run consisted of 128 time points, and two radio frequency excitations were performed prior to image acquisition to achieve steady-state transverse relaxation.  2.2.3  fMRI Pre-Processing  The fMRI data of each individual was pre-processed independently for motion correction, smoothing, and time realignment using Statistical Parametric Mapping (SPM-99) [6]. The time series of functional images were aligned for each slice in order to minimize the signal changes related to small motion of the subject during the acquisition.  Spatial filtering of functional time series was performed by con-  volving each E P I image with a two-dimensional (2D) Gaussian smoothing kernel with full width at half maximum ( F W H M ) of 2.8 mmx2.8 mm. Temporal filtering of functional time series included removal of the linear drifts of the signal with respect to time from each voxel's time-course and low-pass filtering of each voxel's time-course with a one-dimensional (ID) Gaussian filter with F W H M of 6 s. After preprocessing, a voxel-based T-statistical map was generated for each individual by correlating the time course of each voxel in the epochs of interest with an empirically determined hemodynamic response (HDR) obtained from [7].  2.3 2.3.1  Methods Framework of Proposed Method  The overarching goal of the method is to provide sensitive linear discriminant analysis for fMRI activation statistics within a given ROI. This may either be across  46  subjects, where the same ROI is compared across subjects belonging to different groups (e.g. for the ROI, "right cerebellar hemisphere", "disease group" is compared to "control group"), or may represent an ROI in the same group of subjects across different tasks (e.g. for the ROI, "right cerebellar hemisphere", the activations statistics for "internally guided task" are compared to "externally guided task"). For simplicity, in the remainder of this chapter we will simply refer to both these situations as comparing "Group A" to "Group B " . We first assume that T-statistics associated with activation have already been calculated using, e.g. the General Linear Model as utilized in SPM, and that "activated" voxels have been determined by finding those voxels exceeding an arbitrary threshold (e.g. 1.96). Fig. 2.1 summarizes the framework of the proposed method. The 3D moment invariants of the activated-voxel statistics for each subject/condition are then collected into two groups. The invariants across all groups are then dimension-reduced and projected onto a space with maximum discrimination using linear discriminant analysis (LDA). A permutation test is then applied to the resulting feature vectors to estimate the level of significance (p-value) of the null hypothesis that, the activation distributions of the two groups are the same. The implementation details of each module are described in Sections 2.3.2 to 2.3.4.  2.3.2  3D Moment Invariants  The 3D moments of a density function p(x,y,z), (representing functional activation levels in the present context), are defined as: +oo  /  r+oo r+oo  / -oo J— oo  /  x y z p(x,y, p  J—oo  47  q  r  z)dxdydz  (2.1)  T-statistics Group A Group B  , I  1 ,  3D Moment Invariants  LDA  Permutation Test  I  p-value  Figure 2.1: Framework of proposed method. where  n = p+q+r,  is the order of the moment,  (x,y,z)  are the coordinates of each  voxel, and p(x,y,z) are the T-statistics representing the level of activation of each voxel inside the ROI. Like the one dimensional (ID) case where the moments describe the shape of the density function, 3D moments, similarly, provide the general 3D shape information of p(x,y,z). However, potential spatial misalignments in the T-maps would be problematic to the analysis if the 3D moments were naively applied directly, as the misalignments would be inappropriately interpreted as shape changes. Therefore, metrics invariant to rotational, translational, and scaling artifacts are required. In this study, 3D moment invariants were used to ensure that only true shape differences in the activation are captured, independent to the particular co-ordinate system used. Translational invariance can be obtained by using central moments defined as: +oo  /  r+oo  / •oo  J—oo  r+oo  /  (x - x) (y p  -y) {z q  - zf  • p(x,y, z)dxdydz  (2.2)  J—oo  where x, y, and z are the centroid co-ordinates of the density function, calculated 48  as: x —  ,y — mooo  ,z —  .  rnooo  l^-^J  mooo  Scale invariance can be obtained by normalizing the moments as follows:  (2.4) 3  T  1  m,000 3  To obtain rotational invariance, the 3D moments need to be summed in a certain fashion as given in (2.5) to (2.7) and (2.13) to (2.16). Seven invariants were used in this study, three of which are based on 2  n d  order moments derived by Sadjadi and Hall [8]:  J\ = ^200 + "1020 + m 02  (2.5)  0  J2 = m oo"i<02o + m oo"ioo2 + ^020^002 - m\ 2  2  - mf  0l  J3 = m. oomo omoo2 - ^002^110 + 2 m i i m i o i m o i i - m 2 0 ^ i o i 2  2  - m l  1 0  _  0  0  n  (2.6)  ^200^011 ( - ) 2  7  The remaining four higher order invariants used are based on moment tensor contraction [9]. Considering only affine transforms, general tensors become Cartesian tensors, thus covariant and contravariant tensors are equal [10], meaning if ay is a Cartesian tensor, then: ay = a) = a .  (2.8)  ij  By contracting a £ J „ , a ^  , a j g g , and a% £? k  o p  , 3, rd  4, th  5, th  and 6  th  order 3D  moment invariants were derived. Below we detail the algorithm for generating a 3  rd  order 3D moment invariant based on moment tensor contraction: Let M  l  m  n  be a 3  r d  order moment tensor as defined in [11], where I, rn, n =  1, 2, 3, since 3D moment invariants were pursued. We can obtain a 3  49  r d  order 3D  moment invariant, MI3D,  a£J = a  i j k l m n  n  by contracting  =MM ijk  =M  lmn  where N = (order + 1) • + ° : 2  l  m  a? : mn  M  n  l  m  n  = E ^ a *•m  piqiri  3D  = 10 and p , q  = (3 + 1) • ^  d e r  = MI ,  2  t  if  (2.9)  and n can be  calculated as follows:  Pi(l, m, n) - 5(1 - 1) + £ ( m - 1) + <5(n - 1),  (2.10)  qi  (l, m, n) = 6(1 - 2) + 5(m - 2) + 5(n - 2),  (2.11)  n(l, m, n) = 5(1 - 3) + 6{m - 3) + <5(n - 3),  (2.12)  ai is calculated by counting the number of times the same pattern of pi,  rj is  generated by the above algorithm. For example, a\ corresponding to 772300 equals to 1 since in only the case when I — 1, m — 1, and n — 1 can the pattern {3, 0, 0} be generated. The derived higher order moment invariants used in this study are as follows:  #3  =  m  300 +  m  + 3m?  B  A  =  77-4,0 +  030 + ™003 +  + 4m^  ™040 +  #5 =  "7500  31  004 +  m  + 6m^  u  + 6m§  02  311  210 +  + 10m  4  m  310 +  + 4m?  + 4mg  2 2  + ™050 + ™005  + 20m  m  3  ™201 +  3  ™120 +  6  m  l l l  (2.13)  + 3mo2i + 3 ^ 0 1 2  0 2  + 12m|  3  302  +  5 m  4  m  301 +  + 12m?  30  21  6  m  220  + 12m?  12  + 4m?  (2.14)  03  13  410 +  + 10m|  30  5 m  4 0 1 + 10m§  2 0  -I- 3O777I21 + 30m2  12  + 10m^  03  (2.15) +  5777x40 + 2 0 7 7 7 i  +  !0"7032 +  3 1  + 30777,122  !0 023 + m  5  m  + 20m?  014  50  13  +  5777,i  04  +  5?77o4i  •Be = mloo  + m§ + mg + 6m^ 60  06  10  6m§ x +  +  I5m  2  0  420  + 30m + 1 5 m  402  + 20m3  + 15m  + 60m  231  + 90m 2 + 60m 3 + 1 5 m  + 30m + 6 0 m  132  + 60m 3  411  240  141  + 15m^  42  + 20m^  33  30  + 60m3  22  21  + 60m3 + 20m3 12  21  12  + 30m + 6 m 114  + 15mQ + 6moi 24  10  + 6m  03  (2.16)  2  204  5 + 6m  l50  0 5 1  5  Other higher order invariants can be derived using similar procedures to increase the number of features which might enhance discriminative power of the feature vector. However, the chosen set of 3D moment invariants provided adequate shape discriminative capability for the purposes of this study, thus higher order invariants were not used. 2.3.3  Linear Discriminant Analysis (LDA)  Once the 3D moment invariants of the activation distributions of each subject are calculated, the next step is to test whether those feature vectors (corresponding to the spatial distribution of the activation regions) differ between the two experimental groups. To maximize discriminability, the 3D moment invariants are projected onto a ID space, where the between group variance is maximized, and the within group variance is minimized. This projection is achieved using LDA [12]. The cost function used in LDA is the generalized Rayleigh quotient defined as: wSw  .  T  B  ,  .  where w is the vector pointing in the direction of maximum discrimination of the data, and SB and Sw are the between-group and within-group scatter matrices defined as: SB = (mi - m ) ( m i - m )  t  2  51  2  (2.18)  S  w  = ^Ui^xeCiiix  - rrti)  ~ mi){x  (2.19)  where rrii is a vector containing the mean of each 3D moment invariant of group i, a; is a matrix containing the 3D moment invariants of all subjects, and fij denotes group i. The w that maximizes J(w) must satisfy:  SBW = XSyyU!  (2.20)  for some constant A. By transforming J(w) into a generalized eigenvalue problem, the potential singularity problem associated with SB and Sw (due to limited number of samples) is mitigated. The vector w can be found by determining the eigenvector corresponding to the maximum eigenvalue of (2.20).  2.3.4  Permutation Test  To determine the probability of the activation distributions of the two experimental groups being the same, permutation was used.  The permutation test does not  require a priori assumptions about the data distribution, and thus is preferred over the T-test and F-test [13]. The procedures of the permutation test algorithm are as follows [14]:  1. Given m n-dimensional feature vectors for each group, X and Y where m > n, calculate the Mahalanobis distances between each vector in X and the centroid of Y as follows:  (2.21)  where y is the centroid of Y, x j is the j h element of the %th sample feature l  t  vector of X, and Sk is square root of the k  th  covariance matrix of X and Y. 52  diagonal element of the pooled  2. Calculate mean of ^Ma/iaKO and denote this mean as d i . or  g  3. Pool all sample feature vectors, stripped of their group labels. 4. Randomly draw m samples from a random distribution and label these samples X' and label the remaining samples Y'. 5. If the mean dMahal(i) between X' and Y' is > d i or  g  increment counter q by 1.  6. Repeat steps 3 to 6 (N-l) times (ex. N = 5000). 7. The p-value representing the significance of the differences between X and Y is estimated by q / N.  2.3.5  Visualization of Activation Patterns  Since the voxel distributions are 3D, the T-statistics of each voxel can only be represented by color or intensity.  However, plotting the T-statistics in this form  does not provide much insight into the spatial distribution of the activation, since shapes in the 4  th  dimension are hard to visualize. Therefore, a different approach is  taken. Projecting the 3D structural data onto a 2D plane through P C A and using the T-statistics as the 3  r d  dimension, surface plots of the activation distribution were  generated through cubic interpolation. Note that the first two principle components were used in generating these surface plots to ensure that most of the structural information is retained Fig. 2.2.  2.3.6  Comparison of ROI Analysis Methods  Given the T-values of each voxel inside an ROI, metrics for representing the characteristics of the activation within the ROI are required so that comparisons can be 53  Left prefrontal ROI Externally-guided task  3D Moment Invariants f1  I <2  |O |M  I »5 I f$ |  f?  Intensity b a s e d \ methods  Thresholded Mean  Internally-guided task  j  3D Moment Invariants  I fi I f2 I  f3  I H I f5  116 117  Percentage Activated Voxel  I  Figure 2.2: Methods used to compare activation statistics across ROIs. The proposed method captures the rich shape information of the activated regions using 3D moment invariants. Note that irrespective of the orientation of the activation (as demonstrated by the 2 different perspectives shown) the same 3D Moment Invariants will be generated. Conventional methods (right) collapse the spatial information into a single metric such as the mean thresholded T-statistic (upper right) or percentage of activated voxels (lower right), ignoring the texture of the activation, and potentially reducing discriminability across subjects.  54  made between subjects and across groups under different experimental conditions. Conventional methods involve using mean T-values and percentage of activated voxels as metric, whereas we propose incorporating spatial information of the activated voxels. Methods based on these three metrics were compared for detecting the activation differences during externally and internally guided tasks: • T h r e s h o l d e d M e a n : Taking T-values greater than 1.96, compute the mean T-values during externally and internally guided tasks. Calculate the difference between the means of these two conditions, and apply T-test to determine if the differences are significant. • % A c t i v a t e d Voxels: Compute the percentage of voxels with T-values greater than 1.96 during externally and internally guided tasks. Calculate the difference between the percentages of these two conditions, and apply T-test to the differences to determine level of significance. • 3 D M o m e n t Invariants: Proposed method as described in Section 2.3.1. The first two methods are intensity-based, whereas the 3D Moment Invariant method incorporates the extensive spatial information of the activated voxels. In this study, 10 ROIs were examined, consisting of the left and right supplementary motor areas (SMAs), cerebellar hemispheres, primary motor cortices, prefrontal cortices, and caudates. The activation distributions of these ROIs during externally versus internally guided finger sequential movements were compared.  2.4  Results and Discussion  The results generated using the three methods described in Section 2.3.5 are summarized in Table 2.1. As suggested by Fig 2.2, conventional methods, by collapsing the 55  Table 2.1: p-values of Activation Differences During External versus Internal Task Experimental Conditions ROIs Left S M A Right S M A Left C E R Right C E R Left P M C Right P M C Left P F C Right P F C Left C A U Right C A U  Thresholded Mean 0.0466* 0.2955 0.5853 0.7510 0.1496 0.6656 0.8641 0.7140 0.1423 0.6369  Right Hand % Activated Voxels 0.2134 0.0880 0.7093 0.8793 0.5681 0.9110 0.7796 0.0741 0.0469* 0.0757  3D Moments Invariants 0.0352* 0.0796 0.0158* 0.0174* 0.3188 0.1076 0.0334* 0.0328* 0.1014 0.0048*  Thresholded Mean 0.0328* 0.0454* 0.0828 0.1280 0.4301 0.1027 0.2537 0.5955 0.0440* 0.2358  Left Hand % Activated Voxels 0.0796 0.0178* 0.6899 0.3824 0.3243 0.9874 0.1209 0.3289 0.0427* 0.1178  3D Moment Invariants 0.0668 0.0238* 0.1602 0.1316 0.0056* 0.0738 0.0378* 0.3208 0.1600 0.0038*  C E R = cerebellum, P M C = primary motor cortex, P F C = prefrontal cortex, and C A U = caudate. * Statistically significant at a = 0.05.  T-statistics within an ROI into a single value has the advantage of being invariant to the co-ordinate system used and enables values across subjects/conditions to be combined in a straightforward manner. Nevertheless, these advantages are mitigated by the loss of potentially useful spatial information. Unless the spatial information - in effect, the 3D texture - of an ROI is utilized, the changes in T-statistics during externally and internally guided tasks might be very difficult to detect using only intensity information (Fig. 2.2, right). On the other hand, examining the shapes of the activation distributions from different angles reveals that the patterns of activation are substantially different during externally and internally guided tasks. Thus, exploiting this spatial information may provide higher sensitivity to activation differences, provided these differences are not subject-specific. A key result of this study is that the enhanced sensitivity of the 3D moments does not necessarily amplify inter-subject differences to the same extent. We found consistent, significant differences in the spatial distribution of ROI-based activation statistics even when combining across subjects (Table 2.1). This suggests that some 56  spatial aspects of fMRI activation within ROIs are modulated by the task performed, and do not simply reflect inter-subject variability. There is increasing interest in the role of BOLD signal decreases reflecting neuronal de-activations (e.g. [15]). Although not examined in the current study, we note that the proposed method can easily incorporate negative T-values, by concatenating 3D moment invariants associated with the negative T-values into the feature vector for classification. 2.4.1  Supplementary Motor Areas  Both the thresholded mean method and the 3D invariant method detected activation differences between externally and internally-guided tasks in the left SMA when the dominant hand was used. Since the SMA is known to be preferentially invokedduring internally guided tasks [16] [17] [18], detecting activation differences in the left SMA is expected, based on prior neuroscience knowledge, when the task involves the right hand. During the left handed task, all three methods detected activation differences in the right SMA. This again matches expectation since the SMA is more involved during internally guided tasks, and the left hand is used. However, the thresholded mean T-statistic also detected significant activation differences in the left SMA. This might have arisen from connections between the left and right SMA, but requires connectivity analysis to verify.  2.4.2  Cerebellum  Only the proposed method was able to detect significant activation differences in the cerebellum during the right handed task. This activation difference might be  57  associated with handedness as seen in past studies [19].  2.4.3  Primary Motor Cotex  Only the 3D moment invariant method detected differences using the non-dominant hand. As recent evidence has suggested distinct mechanisms are utilized with the dominant and non-dominant hand [20], [21], we suggest that the different activation between tasks with the non-dominant hand is expected.  2.4.4  Prefrontal Cortex  Again, only the proposed method was able to detect significant activation differences in both prefrontal cortices during the right handed task and in the left prefrontal cortex during the left handed task. Fig. 2.2 suggests why only the proposed method detected activation differences (i.e.  the mean T-value and percentage of  activated voxels differ by a minuscule amount during the two tasks, whereas substantial changes in the shapes of the activation distributions are shown in the Fig. 2.2).  Prefrontal cortex has been previously shown to be differentially activated  during externally guided and internally guided tasks [22].  2.4.5  Caudate  Both intensity-based methods detected activation differences in the left caudate during the left handed task. In addition, percentage activated voxel method detected activation differences in the left caudate during the right handed task. 3D moment invariant method, on the other hand, detected no activation differences in the left caudate, but detected activation differences in the right caudate in both tasks. Since the right caudate is associated with timing information [23] and the load placed on  58  the right caudate to process the timing information during the two tasks are quite different, the detected activation differences is expected.  2.5  Conclusion  We have demonstrated that certain spatial aspects of activation statistics appear to be relatively well conserved across subjects yet are significantly modulated by task performance. Incorporating spatial information of activation statistics within an ROI by utilizing 3D moment invariants appears to complement traditional activation based methods. A direct extension of the methods is to look at the spatial characteristics of combinations of ROIs that maximally discriminate between groups, an approach currently being pursued.  59  Bibliography [1] J. L. Lancaster, L. H. Rainey, J. L. Summerlin, C. S. Freitas, P. T. Fox, A. C. Evans, A. W. Toga, and J. C. Mazziotta. Automated labeling of the human brain: A preliminary report on the development and evaluation of a forwardtransform method. Human brain mapping, 5(4):238-242, 1997. [2] G. Allen and E. Courchesne. Differential effects of developmental cerebellar abnormality on cognitive and motor functions in the cerebellum: an fmri study of autism. The American Journal of Psychiatry, 160(2):262-273, Feb 2003. [3] H. J. Aizenstein, K . A. Clark, M . A. Butters, J. Cochran, V . A. Stenger, C. C. Meltzer, C. F. Reynolds, and C. S. Carter. The bold hemodynamic response in healthy aging. Journal of cognitive neuroscience, 16(5):786-793, 2004. [4] P. Bogorodzki, J. Rogowska, and D. A. Yurgelun-Todd. Structural group classification technique based on regional fmri bold responses. IEEE  Transactions  on Medical Imaging, 24(3):389-398, Mar 2005. [5] J. F. Mangin, F. Poupon, E. Duchesnay, D. Riviere, A. Cachia, D. L. Collins, A. C. Evans, and J. Regis. Brain morphometry using 3d moment invariants. Medical image analysis, 8(3): 187-196, Sep 2004.  60  [6] http://www.fil.ion.ucl.ac.uk/sprn/. [7] S.A. Huettel and G. McCarthy. Evidence for a refractory period in the hemodynamic response to visual stimuli as measured by MRI. Neurolmage, 11(5):547553, 2000. [8] F. A. Sadjadi and E. L. Hall. Three-dimensional moment invariants.  IEEE  Transactions on Pattern Analysis and Machine Intelligence, 2:127-136, 1980.  [9] T. H. Reiss. Features invariant to linear transformations in 2d and 3d. pages 493-496. 11th IAPR International Conference onPattern Recognition, The Hague, Netherlands, 1992. [10] D. Cyganski and J. Orr. Object recognition and orientation determination by tensor methods. Advances in Computer Vision and Image Processing, 3:101-  144, 1988. [11] D. Cyganski and J. Orr. Applications of tensor theory to object recognition and orientation determination. IEEE  Transactions on Pattern Analysis and  Machine Intelligence, 7(6):662-673, 1985.  [12] R. O. Duda, P. E. Hart, and D. G. Stork. Pattern Classification.  Wiley-  Interscience, 2000. [13] J. Ludbrook and H. Dudley. Why permutation tests are superior to t and f tests in biomedical research. The American Statistician, 52(2): 127-132, 1998. [14] S. Vetsa, M . Styner, S. Pizer, J. Lieberman, and G. Gerig. Caudate shape discrimination in schizophrenia using template-free non-parametric tests, pages 661-669. Medical Image Computing and Computer-Assisted Intervention-  61  MICCAI 2003 6th International Conference, Montreal, Canada, Nov 15-18, 2003. [15] J. Pfeuffer, H. Merkle, M . Beyerlein, T. Steudel, and N . K . Logothetis. Anatomical and functional mr imaging in the macaque monkey using a vertical largebore 7 tesla setup. Magnetic resonance imaging, 22(10):1343-1359, Dec 2004. [16] J. M . Orgogozo and B. Larsen. Activation of the supplementary motor area during voluntary movement in man suggests it works as a suprarhotor area. Science (New York, N.Y.), 206(4420):847-850, Nov 16, 1979.  [17] P. E. Roland, B. Larsen, N . A. Lassen, and E. Skinhoj. Supplementary motor area and other cortical areas in organization of voluntary movements in man. Journal of neurophysiology, 43(1):118-136, Jan 1980.  [18] G. Goldberg. Supplementary motor area structure and function: review and hypotheses. Behavioral and brain sciences, 8(4):567-615, 1985.  [19] A. Solodkin, P. Hlustik, D. C. Noll, and S. L. Small. Lateralization of motor circuits and handedness during finger movements. European journal of neurology : the official journal of the European Federation of Neurological Societies,  8(5):425-434, Sep 2001. [20] R L Sainburg and D. Kalakanis. Differences in Control of Limb Dynamics During Dominant and Nondominant Arm Reaching. Journal of Neurophysiology, 83(5):2661-2675, 2000. [21] L. B. Bagesteiro and R. L. Sainburg. Handedness: dominant arm advantages in control of limb dynamics. Journal of neurophysiology, 88(5):2408-2421, Nov 2002. 62  [22] I. H. Jenkins, M. Jahanshahi, M. Jueptner, R. E. Passingham, and D. J. Brooks. Self-initiated versus externally triggered movements, ii. the effect of movement predictability on regional cerebral blood flow. Brain : a journal of neurology, 123(6):1216-1228, Jun 2000. [23] R. Jech, P. Dusek, J. Wackermann, and J. Vymazal.  Cumulative blood  oxygenation-level-dependent signal changes support the 'time accumulator' hypothesis. Neuroreport, 16(13):1467-1471, Sep 8, 2005.  63  Chapter 3  Feature Selection with Modified LDA 3.1  Introduction  One main application of functional magnetic resonance imaging (fMRI) is to make statistical inferences about regions of the brain responsible for different behavioral tasks. Typically, this entails assigning a statistical value to each voxel related to its probability of being activated during the performance of the task with all voxel statistics assembled into a statistical parametric map.  The individually-derived  fMRI statistical maps are then either compared across subject groups or across different tasks to draw statistical inferences about group activations. Due to intersubject variability in brain shape, brain size, and orientation of the subjects in the M R scanner, combining fMRI statistics across subjects to make group inferences is A version of this chapter is under revision. B. Ng, R. Abugharbieh, X. Huang, and M . J. McKeown, "Characterization of fMRI Activation Maps Using Invariant 3D Moment Descriptors," I E E E Transactions on Medical Imaging, invited paper (under revision). 2  64  not trivial. To this end, two broad approaches have been traditionally employed. The most common approach involves warping the brain of each individual subject to a common exemplar template [1]. However, spatial normalization, typically followed by spatial smoothing, may inappropriately pool responses from functionally dissimilar regions [2], especially for small, subcortical structures, thus degrading important spatial information. The alternative approach to whole brain warping involves drawing the regions of interest (ROIs) for each individual subject, and examining the statistical properties of regional activations across subjects within those ROIs. Castanon et al. demonstrated that this latter approach offers finer localization and increased sensitivity to task-related effects [2], thus this subject-specific ROI-based approach is adopted in this thesis. Defining a metric to rigorously characterize and compare ROI activation maps across groups of subjects without merely amplifying inter-subject variability is also very challenging. Mean activation statistics, percentage of activated voxels, and mean voxel time course averaged over an ROf have all been previously used [3].  For example, Allen and Courchesne used the percentage of activated voxels  and mean percent signal changes to study the activation differences between subjects with and without autism [4]. Aizenstin et al. used the mean intensity time courses of the significantly activated voxels and the percentage of voxels with negative T-values to compare activation distributions of young and elderly subjects [5]. Similarly, Bogorodzki et al. used features extracted from mean intensity time curves to discriminate between different groups of marijuana smokers [6]. A l l these group discrimination techniques are essentially based on comparing the magnitudes of the activation statistics with the underlying assumption that the spatial distribution of voxel activation statistics within an ROI is not meaningful. These features will  65  thus hereafter be referred to as 'magnitude-based' features in the remainder of this chapter. Characterizing the activation statistics within an ROI using magnitude-based features provides a simple means for comparing activation maps across different subjects, since these methods are insensitive to variability in brain shape, brain size, subjects' orientation inside the MR scanner, and location of the activated voxels. However, we have previously demonstrated that exploiting the spatial information encoded by the location of the activated voxels within an ROI can significantly increase discriminative power [7]. We note an interesting parallel between magnitudebased features, which ignore the spatial distribution of activation statistics within an ROI, and neural coding methods in the brain. At the neuronal level, the mean firing rate (known as "rate coding") is only one of the many possible ways in which a particular brain response might be encoded. Other models, such as "population coding", suggest that the brain response could be inferred from the overall pattern of activity in a set of (spatially-disparate) neurons [8], where this population coding may also be present at the macroscopic scale of neuroimaging [9]. Similarly, instead of only looking at mean changes in the blood-oxygenation level dependent (BOLD) signal averaged over an ROI, it may be beneficial to also explore the spatial distribution changes of BOLD signals during different task-performances. This is because the spatial distribution of activation statistics within an ROf may actually be altered by task performance even when no change in mean activation statistics or percentage of activated voxels is present as illustrated in Fig. 3.1. In this thesis, we propose a novel method for characterizing the spatial distribution of fMRI activation statistics within an ROI using three dimensional moment invariants (3DMIs). The invariance properties of the proposed descriptors are critical in ensuring that  66  Condition B  Condition A  Figure 3.1: Simplified synthetic example of ROI activation maps illustrating the underlying idea of the proposed method. Assuming 100 voxels are within the ROI and only 3 voxels are active in both conditions. Mean activation statistics = 6/100 and percentage of activation = 3% are the same for both conditions A and B. However, these two tasks are spatially encoded quite differently. This difference in the spatial distribution of activation statistics could be characterized using 3D moment invariants as proposed in this thesis to better discriminate changes in activation patterns between the two conditions.  67  the extracted spatial features are independent of the subject's orientation in the M R scanner. Also, using spatially invariant features accounts for some of the intersubject variability such as differences in brain sizes. A recent study by Mangin et al. demonstrated the discriminative power of 3DMIs on structural MR1 data [10]. In this thesis, we extend the use of such invariant descriptors to functional M R I data. Our proposed method is applied to real fMRI data recorded from twenty one healthy subjects performing externally and internally-guided motor tasks. The actual motor task in both cases is the same (see Section 3.3), but the brain processes involved are different depending on whether the subject is mimicking a visual cue or recalling the movement from memory. This relatively minor but nonetheless crucial distinction makes such a task ideal for testing the sensitivity of different activation discriminant analysis methods.  Using 3DMIs, we demonstrate signifi-  cant task-related activation differences in the left lentiform nucleus, right cerebellar hemisphere, left supplementary motor areas (SMA), left primary motor cortex ( M l ) , and left prefrontal cortex (PFC). In contrast, percentage of activated voxels did not detect any significant tasked-related activation differences, and thresholded mean activation statistics only detected activation differences in the SMAs.  These re-  sults suggest that different tasks might be encoded by the 3D 'spatial texture' of the activation statistics within an ROI in addition to the magnitude and extent of activation.  68  3.2 3.2.1  Materials Study Subjects and Experimental Conditions  All research was approved by the appropriate Institutional and Ethics Review Boards. After obtaining informed consent, twenty one healthy subjects were asked to perform two motor tasks in 20 s blocks with their right hand while in the fMRI scanner: • In the first task, subjects were externally guided (EG), following a finger tapping sequence displayed as a movie of a hand on a screen. • In the second task, subjects were internally guided (IG), performing a finger tapping sequence they had seen on a screen earlier during the first task, but recalling the sequence from memory.  3.2.2  Data Acquisition  Images were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangern, Germany) with a birdcage-type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The participant's head was positioned along the canthomeatal line. Foam padding was used to limit head motion within the coil. High-resolution Tl-weighted anatomical images (3D S P G R , TR=14 ms, TE=7700 ms, flip a n g l e = 2 5 ° , voxel dimensions l . O x l . O x l . O mm, 176x256 voxels, 160 slices) were acquired for co-registration with the functional images. Co-planar T2*-weighted functional images were acquired at each time point using a gradient echoplanar sequence (TR=3000 ms, T E = 3 0 ms, flip a n g l e = 8 0 ° , N E X = 1 , voxel dimensions 3.0x3.0x3.0 mm, imaging matrix 64x64 voxels, 49 slices). Each functional run consisted of 128 time points, and two radio frequency excitations were performed prior to image acquisition to achieve steady-state transverse relaxation.  69  3.2.3  fMRI Pre-Processing and Statistical M a p Generation  The fMRI data was preprocessed for each subject independently for slice timing and motion correction using the standard parametric mapping software ( S P M 2) [11].  Spatial filtering of functional time series was performed by convolving each  E P I image with a two-dimensional Gaussian smoothing kernel with full width at half maximum ( F W H M ) of 2.8 mmx2.8 mm. Temporal filtering of functional time series included removal of the linear drifts of the signal with respect to time from each voxel's time course, and low-pass filtering of each voxel's time course with a one dimensional (ID) Gaussian filter with a F W H M of 6 s. No spatial normalization (brain warping) was performed. After preprocessing, each time point in all runs was coded as a particular type of task (such as rest, right hand E G finger sequence movements (FSM), etc.). To emphasize the proposed spatial analysis, we chose to use the simplest method in generating the statistical parametric maps: a simple T-test applied to each voxel based on the B O L D signal changes in all runs between two tasks: 1) right hand E G F S M versus rest and 2) right hand IG F S M versus rest. The resulting voxel activation statistics were then mapped onto the anatomical ROIs, which were manually drawn by a trained observer according to standard atlases [12].  Regions exam-  ined included bilateral caudate, lentiform nucleus, thalamus, S M A , M l , P F C , and cerebellar hemisphere. More details about the ROI specification procedures can be found in [13].  70  3.3  Methods  3.3.1  Proposed Framework for Spatial fMRI Analysis  The overarching goal of the proposed method is to provide sensitive discriminant analysis of fMRI activation statistics within an ROf. This may either be across subjects, where the same ROI is compared across subjects belonging to different groups (e.g.  for the 'right cerebellar hemisphere' ROI, all subjects in a 'disease  group' are compared to all subjects in a 'control group'), or may be within the same group of subjects across different tasks (e.g.  for the 'right cerebellar hemisphere'  ROI, the activations statistics for the IG task are compared to that of the E G task for all subjects). For simplicity, in the remainder of this chapter, we will simply refer to both situations as comparing 'Group A ' to 'Group B'. Fig.  3.2 summarizes the framework of the proposed method.  Assuming  T-statistics Group A Group B Feature Extraction  Winsorization  LDA/Feature Selection Permutation Test j p-value Figure 3.2: Framework of the proposed method.  that activation statistics have been calculated using e.g. the General Linear Model  71  (GLM) as utilized in SPM, features characterizing the activation statistics within an ROI are first extracted for each subject and collected into two groups. To enhance robustness to potential outliers, a winsorization procedure is performed to generate estimates of the feature values. For the case of single ROI feature analysis using traditional magnitude-based features, linear discriminant analysis (LDA) is not necessary since in a ID feature space, there is no alternative axis orientation for increasing discriminability. Instead, a permutation test is directly applied to the single ROI features to estimate the probability (p-value) of the null hypothesis that the activation maps of the two groups are drawn from the same distribution. For the case of multiple ROI feature analysis using the proposed spatial features, LDA is first applied to all subjects' ROI features to determine the relative weightings of these features in the maximal discriminant ID space. These weightings are then used to sort the features in decreasing level of discriminability. A modified LDA procedure that incorporates leave-one-out cross-validation is then performed to choose the number of features to include in the analysis without over-fitting. A permutation test is then applied to the resulting feature vectors to estimate the probability that the activation maps are from the same distribution. 3.3.2  ROI Spatial Feature Extraction: 3D Moment Invariants  Two types of ROI features are contrasted in this thesis. One is based on activation statistics magnitude, whereas the other is based on the spatial distribution (or '3D spatial texture') of activation statistics within an ROI. More specifically, we have chosen to compare two magnitude-based features, thresholded mean activation statistics and percentage of activated voxels, which are widely-used in ROI-based fMRI analysis [3], against the proposed spatial features, namely 3DMIs. The thresh-  72  olded mean activation statistics are calculated by averaging the voxel activation statistics (or T-statistics as used in this thesis) greater than a particular threshold, T , over an ROI. The percentage of activated voxels is defined as the number of voxels with T-values above r divided by the total number of voxels within the ROI. We have tested a range of r from 1 to 4 at 0.5 increments. The ROIs detected with significant activation differences were consistent for this range of thresholds, thus only results generated with r = 1.96, which corresponds to an uncorrected p-value of 0.05, are included for comparisons with the proposed 3DMI-based method. The 3D spatial features we employed (3DMIs) are calculated based on 3D moments defined as: +oo  r+oo r+oo / / x y z p{x,y,z)dxdydz, / -oo J—oo J—oo p  q  r  (3.1)  where n — p+q+r is the order of the 3D moment, (x,y,z) are the coordinates of a voxel, and p(x,y,z) is the T-value of a voxel located at (x,y,z) within an ROI. We have restricted our analysis to positive T-values, since there is little consensus on the interpretation of negative T-values [14]. As for the choice of T-threshold, we have decided to include all voxels with T-values greater than zero. One might argue that using a T-threshold of zero will include many non-active voxels that are falsely deemed active (false positives) into the analysis. However, false positives will not likely be located at the same position within an ROI across subjects. Hence, inclusion of these voxels will only reduce discriminability of our spatial features. On the other hand, truly active voxels that are corrupted by noise may reside in similar location within an ROI across subjects. If a high T-threshold is set, some of these active voxels will be excluded from the analysis. Therefore, we argue that using a T-threshold of zero provides a suitable, conservative estimate of activation differences, while retaining most of the original information of the ROIs, and is thus  73  appropriate for our spatial activation statistics analysis. Naively using 3D moments in discriminating activation statistics is problematic, since potential spatial misalignments of the T-maps would be inappropriately interpreted as spatial changes. Therefore, metrics invariant to similarity transformations (rotation, translation, and scaling) are needed.  Furthermore, by using  spatially invariant features, inter-subject variability in brain size or slight translational/rotational shifts of the activated region would be accounted for. In this study, 3DMIs are used to ensure that only true spatial differences in activation statistics are captured in the fMRI data analysis, i.e. independent of the subjects' orientation in the M R scanner and brain size. Translational invariance of 3D moments can be obtained by using central moments defined as: +oo  /  r+oo  / -oo  r+oo  /.  (x - x) (y  -y) (z  p  - z)  q  r  • p(x,y,z)dxdydz  (3.2)  J—oo J—oo  where x, y, and z are the centroid co-ordinates of p(x,y,z), calculated as:  _ _ mum x —  ?7i<ooo  _ _ mpio ,  y —  mooo  _ _ mpoi ,  z —  mooo  .  (."J.o)  Scale invariance can be obtained by normalizing the moments as follows:  P<ir ~  J  P+g+r m  !  1  '  (3-4)  000  Examining (3.2), if we consider the distance from the centroid (e.g.  (x — x)) as a  weighting of the T-values p(x,y,z), the T-values of voxels closer to the ROI boundary will be more highly weighted. However, some of the outer voxels may not belong to the ROI, but are included due to manual segmentation errors. Also, T-values of voxels near the activation centroid are less prone to noise as higher signal strength is expected at those voxel locations. Thus, to mitigate manual ROI segmentation 74  errors and to stress voxels near the activation centroid (i.e. to ensure errors in Tvalues of voxels near the ROI boundary are not given higher weighting), the following remapping is performed: 1  (x-x) — ), 2  -.exp{  y = —^=exp(-&^-),  (3.5)  a  JCTy  (T \JZ7T y  1  (z-z)\  (  9 = —fK= P\  Z  WZ  eX  a \/2-K  )>  2o  z  z  where (x — x), (y — y), and (z — z) in (3.2) are substituted by x , y , and z , and a , g  g  g  x  a , and o are set as 2 times the size of the ROI along each direction to account for y  z  inter-subject variability in brain sizes. Values lower than 2 can be used to further emphasize the voxels near the centroid and vice versa, but we empirically found that a value of 2 deweights the outer voxels without totally neglecting them. Also, since this remapping bounds x , y , and z between 0 and 1, which accounts for interg  g  g  subject variability in brain sizes, the traditional way for obtaining scale invariance (3.4) is no longer required. We denote the centralized 3D moments with (x — x), (y ~ y)y  a n  d (z — z) in (3.2) replaced by x , y , and z aS T]pq . g  g  g  r  To obtain rotational invariance, the 3D moments need to be combined in specific ways. Sixteen 3DMIs are examined in this study. Three of which are based on 2  n d  order moments derived by Sadjadi and Hall [15]: J\  J2 •^3 =  =  ??200^020  772007702077002 -  =  ^200 +  + ?7200?7002 +  7700277110 +  ?702O  (3.6)  + ?7002  77020*7002  -  2771107710177011  Vwi - VllO -  -  77o2o7??oi -  77011  7720077011  ( - ) 3  7  (3.8)  Another eleven 3DMIs ( - ^ 2 2 ! - ^ 2 2 2 ' ^ 3 3 ! ^ i i > ^ 2 3 3 > ^ i 2 3 i ^ rived by Lo and Don based on complex moments are also used [16]. In addition, we 75  derived two other higher order 3DMIs (.63 and B4) using moment tensor contraction similar to the approach described in [17] as described next. Considering only rigid transformations, general tensors become Cartesian tensors, thus covariant and contravariant tensors are equal [18], meaning if ciij is a Cartesian tensor, then: ay = oj = a .  (3.9)  ij  Thus, by contracting a]^  and a^nop, 3  n  and 4  r d  order 3DMIs can be derived as  th  follows:  For the 3  r d  order case, let M\  be a 3  order moment tensor as defined in [19],  rd  mn  where each of the subscripts: I, m, and n, goes from 1 to 3 since 3DMIs are pursued. By contracting , we obtain a 3  a^  n  = a  i j k l m n  = M  i j k  M  order 3DMI, B : 3  l m n  where JV = (order + 1) • + ° ^ 2  r d  d e r  = M  l  m  n  M  l  m  n  = E ^ a *•m  = (3 + 1) . 2±3  =  1  0  = B,  2 p i q i r i  and p  u  (3.10)  3  q , and TJ can be {  calculated as follows:  Pi(l,m,n)  = 5(1 - 1) + 5(m - 1) + 6(n - 1)  (3.11)  qi(l,m,n)  = 6(1 - 2) + 6(m - 2) + 5(n - 2)  (3.12)  n(l, m, n) = 5(1 - 3) + 5(m - 3) + 5(n - 3)  (3.13)  where is the 5(-) Kroneker's delta function. Q i is calculated by counting the number of times the same pattern of pi, qi, T-J is generated by the above algorithm. For example, a.\ corresponding to  771300  equals  to 1 since in only the case when I — 1, rn = 1, and n — 1 can the pattern {3, 0, 0}  76  be generated. Based on the above tensor contraction procedures, the resulting 3 order 3DMI is defined as below: #3 = ?7300 + ??030 + 77003 +  3^7210 + ?20i + 7i20 + 7 i i i 3r  3r  6r  + 3T7?2 + 3r?o2i + 3r?oi2 0  (3-14)  BA = T7400 + 77040 + 7?004 + 7?310 + 7301 + ''7220 + 1277211 + 677202 + 4T7? 4  4T  6  30  (3.15) + 12i7? + 1277^ + 477? + 477^ + 677^22 + 21  03  In fact, other higher order invariants can be derived using similar procedures. However, to ease interpretation, we have decided to restrict our analysis to 3DMIs that are combinations of 2 , 3 , and 4 order 3D moments, which corresponds to nd  rd  th  combinations of spatial variance, skewness, and kurtosis. For example, J\ can be interpreted as the sum of spatial variance along the x, y, and z direction, and £?3 can be interpreted as a non-linear combination of skewness and other spatial descriptive terms. Note also that higher order moments are less robust to noise [20], thus are not used in this study. 3.3.3  Feature Winsorization  Many linear procedures such as LDA are prone to outliers. To enhance robustness, we first winsorized [21] all the feature values for each ROI separately prior to subsequent processing or analysis. The winsorization procedure is as follows: Let F be a vector with each element, Fi, corresponding to a single ROI feature value (e.g. J\) of the i subject in e.g. group A. Fi of all subjects from both groups are pooled th  into F. To winsorize F, wefirstdetermine which subject(s) are outliers using (3.16): lout  where  J t ou  = find{\  F- F  med  is the set of outlier subjects,  F \ mea  |> N  std  • a),  (3.16)  is the median of F, and N d is the st  number of standard deviations from F ^ beyond which Fi is considered an outme  77  lier. Assuming F, is normally-distributed (since LDA is derived assuming the input features are Gaussian), N  is set to 3 to retain 99.7% of the possible F; values as  std  non-outliers, a is a robust estimate of the standard deviation of F calculated as: a w A • median(\F  - F  |),  med  (3-17)  where A converts the median of the absolute differences (MAD) of F to standard deviation and is set to 1.483 assuming F{ is normally-distributed [21]. For only the outlier Fj's, we winsorize each of these values using (3.18): Ff  n  = F  med  + sgn(Fi - F ) med  •N  std  -a,ie I . mt  (3.18)  This winsorization procedure is applied to each ROI and each feature separately. We note that for the data examined in this study, typically only 1 or 2 subjects were detected as outliers for each feature. Also, we have chosen to winsorize as opposed to throwing away the outliers for two main reasons. First, the number of subjects in most fMRI studies is typically quite low (i.e. like 8 or 10), which renders it impractical to eliminate a given subject's data just because it is an outlier for a particular (single) ROI or feature. However, one might argue that we can always just remove a subject for a particular ROI or feature (i.e. retaining that subject when analyzing the other ROIs and features). The problem is that changing the sample size across ROIs will introduce variability in statistical power, thus fair comparisons across ROfs cannot be made (i.e. not detecting significant differences in activation statistics with a certain ROI could be attributed to reduced number of subjects). Therefore, to avoid complications arising from changes in sample size, we have chosen to winsorize instead of eliminating outliers from further analysis.  78  3.3.4  Feature Selection: Modified L D A  A major issue in fMRI study is the low number of subjects available for making statistical inferences.  If we naively include all described 3DMIs in the analysis,  overfitting is likely to occur. Therefore, we developed a modified L D A procedure that incorporates leave-one-out cross-validation for choosing which 3DMIs to include and in validating our results. Before feature selection can be performed, the 3DMIs must first be sorted based on their overall relative discriminability across ROIs. For every subject and for a particular ROI, we concatenate all the 3DMIs into a single vector, and apply L D A [22] to the resulting feature vectors of the two groups. To mitigate potential singularity problems associated with the within-group scatter matrix, Sw> the projection matrix, W, is determined by solving the following generalized eigenvalue problem: SW B  = XS W, W  (3.19)  where SB is the between-group scatter matrix and A is a diagonal matrix with its diagonal elements corresponding to the eigenvalues of the columns of W. The projection vector, w, that points in the direction of maximal discriminability is then determined by finding the column of W associated with the maximum eigenvalue. Since | w | can be interpreted as the relative weightings of the 3DMIs for which the two groups of feature values are maximally discriminated, the 3DMls can be ranked in descending level of discriminability based on the magnitude of | w \. In this study, 16 3DMIs are examined, so we first rank the 3DMIs from 1 to 16 for each ROI (16 corresponding to the largest | w |). We then sum up the ranks of each 3DMI across ROIs and sort the 3DMIs in descending order of their rank sums. The reason for ordering the 3DMIs based on the overall discriminability across ROIs is to ensure 79  that the same features are used to analyze each ROIs. Else, fair comparisons across ROIs cannot be made. Once the 3DMIs are sorted in decreasing level of discriminability, we employ a modified L D A procedure that incorporates leave-one-out cross-validation to determine the number of 3DMIs, K, to be used in the analysis. The procedure is as follows: 1. Let k be the number of features (ordered in decreasing level of discriminability) to be included in each subject's feature vector, and initialize k to 2. 2. Apply L D A to the two groups of feature vectors, but with the i  th  subject from  both groups excluded, to determine w. 3. Project the i  th  features as F  A  subject's feature vectors onto w, and denote the projected and  Ff.  4. Repeat steps 2 and 3 for every other subject, vertically concatenate the projected features (i.e. Ff- and Ff),  and denote the resulting column vectors as  Ff- and F f . 5. Calculate the Mahalanobis distance between Ff  = Vd•S-  d ahal(K) M  d = mean(F  A  E = Cov(F  A  1  - F ), B  - F ), B  and Ff  • d,  using (3.20):  (3.20)  (3.21) (3.22)  6. Repeat 2 to 5 for every ROI, and denote the sum of all ROIs' d^ahai as d^} . ahal  7. Repeat steps 2 to 6 for k — 3 to 16.  80  8. The k corresponding to the maximum  * * s  dj$a/iai>  n e  n u m  b e r of features  that best discriminates the two groups of activation maps across different ROIs in general. We note that since in step 2, the i  th  subject's information is not used in generating  w, the projected feature values calculated in step 3 are "test data" and not "training data", thus provides validation to our results. A plot of  d^l  ahal  versus the number  of features is provided in the Section 3.4. Also note that the K most discriminant 3DMIs are first projected onto the maximal discriminant ID space using the above modified L D A procedure prior to applying permutation test.  3.3.5  Permutation Test  To determine the probability that the features extracted from the ROI activation maps of the two experimental groups are drawn from the same distribution, we decided to use the permutation test.  Permutation test does not require a priori  assumptions about the data distribution, and is thus preferred over T-test and F test [23]. The permutation test procedure is as follows [24]: 1. Given m n-dimensional feature vectors (assuming row vectors) for each group, vertically concatenate the m feature vectors of the two groups into matrices, X and Y. In the present study, m = 21 (i.e.  21 subjects) and n — 1 (i.e.  each subject's 3DMI feature vector has been projected onto a ID space during L D A , and mean activation statistics magnitude and percentage of activation are ID features). 2. Calculate the Mahalanobis distances between X and Y using (3.20), and denote it as dorig.  81  3. Pool all sample feature vectors, stripped of their group labels. 4. Randomly draw m samples from the pool of feature vectors, label these samples X', and label the remaining samples 5. If  dMahal  between X' and Y' is > d i ,  6. Repeat steps 3 to 6 (N-l)  or g  Y\ increment a counter q by 1.  times, e.g. let N = 100,000.  7. The p-value representing the probability that X and Y are drawn from the same distribution is estimated with the value q/N.  3.3.6  Comparison of R O I Features  As mentioned in Section 3.3.2, two types of features are examined in this thesis: magnitude-based features including thresholded mean activation statistics and percentage of activated voxels, and spatial texture-based features, 3DMIs. We combined the described 3DMIs using the proposed method, and compared the results with that generated using thresholded mean activation statistics and percentage of activated voxels as single features as well as with these two magnitude-based features combined using the proposed analysis framework. We note that one might argue that any enhanced discriminability with the proposed method can easily be due to more features being used. However, in past ROI-based fMRI studies, only a limited number of magnitude-based features were considered due to inter-subject variability in brain shape, brain size, and subject's orientation in the M R scanner. Hence, we are extending the prior set of features with invariant spatial descriptors to make multiple feature analysis feasible.  82  3.4  Results  In this study, we are interested in comparing magnitude-based features (i.e thresholded mean activation statistics and percentage of activated voxels) against the proposed spatial features, 3DMIs. To determine whether the differences in activation statistics during E G versus IG are significant, permutation test was applied to the described features with the results summarized in Table 3.1.  Table 3.1: p-values of Activation Differnces During External versus Internal Tasks based on Positive T-statistics , , ROI  Thresholded Mean  Left L E N Left C E R Left SMA Left T H A Left M l Left P F C Left C A U Right L E N Right C E R Right SMA Right T H A Right M l Right P F C Right C A U  0.1565 0.3143 0.0036* 0.2822 0.1709 0.5254 0.4146 0.4559 0.5079 0.0416* 0.2674 0.8176 0.6828 0.4587  Percentage of Activated Voxels 0.5576 0.7178 0.1461 0.2902 0.5256 0.6588 0.3759 0.9640 0.2656 0.0709 0.0942 0.7541 0.3062 0.9291  Proposed 3DMIs 0.0383* 0.4412 0.0184* 0.1144 0.0409* 0.0469* 0.8351 0.3330 0.0252* 0.9895 0.3465 0.3794 0.8906 0.2097  L E N = lentiform nucleus, C E R = cerbellar hemisphere, S M A = supplementary motor area, T H A = Thalamus, M l = primay motor coretex, P F C = prefrontal cortex, and C A U = caudate. *statistically significant at a = 0.05.  The percentage of activated voxels did not detect any.significant activation differences in any of the ROIs, and the thresholded mean only detected significant differences in the SMAs. In contrast, by first arranging the 3DMIs in decreasing order of discriminability lf , lf , 133  I ,J2,Ji 2  2  n3  if  3  3  3  , iff , 3  I l , / | , JzJn, 2  2  3  B , I% ,7|| , 3  333  3  I$ , B , 22  4  and using the proposed modified L D A procedure to choose the optimal  83  number of features (K = 5 as shown in Fig. 3.3), the proposed 3DMI method was able to detect significant task-related differences across subjects in the left lentiform nucleus, right cerebellar hemisphere, left SMA, left M l , and left PFC. 4.500 4.000 3.500 |  3.000  I 2.500 .2 |  2.000  a  o 1.500 1.000 0.500 0.000 -I  1  1 1  2  1 3  1 4  1  1 5  6  1  1  1  7 B 9 10 Number of Features  1  1 11  1 12  1 13  1  1 14  15  16  Figure 3.3: Spatial feature selection (see description in Section 3.3.4). With the 3DMIs ordered based on their relative discriminability, 5 appears to be the optimal number of features for the dataset examined in this study. When more than five 3DMIs are used, discriminability is reduced due to overfitting.  3.5  Discussion  Our results suggest that the proposed 3DMl-based method is much more sensitive in detecting task-specific group changes in fMRI activation statistics within an ROI than using magnitude-based features. The ROIs detected with task-specific differences are consistent with previously-described neurobiological knowledge, as briefly 84  summarized below.  3.5.1  Neurobiological Relevance of Task-specific ROI Differences  The lentiform nucleus is known to be involved in diseases such as Parkinson's disease, where subjects can more easily perform E G tasks [25],[26],[27]. This notion is further supported by a recent study which reported increased activity in the basal gangliathalamo-cortical loop during IG tasks that was not present during E G tasks [28]. Thus, detecting significant differences in the lentiform nucleus is consistent with prior imaging findings. Significant differences in activation statistics were also detected in the right cerebellar hemisphere. This difference might be associated with handedness as seen in past studies [29]. Previous research has also indicated that sequential (as opposed to single) tasks often recruit the contralateral cerebellar hemisphere [29], whereas the typical motor pathway involves the ipsilateral cerebellar hemisphere [30]. It has been postulated that this increased activity in the contralateral cerebellar hemisphere may be related to the encoding of complex tasks [31]. Detecting significant differences in the left S M A can be explained by the fact that this neural region is known to be preferentially invoked during internally guided tasks [32], [33], [34]. This difference attributes to the increased contralateral activity for right-hand tasks. Significant differences between E G and f G tasks were also detected in the left M l and left P F C . The failure of magnitude-based methods in detecting significant differences in the left M l could be explained by the fact that the M l is the final pathway for motor movements and the actual movements performed during the E G and IG conditions were similar, thus requiring about the same level of activation.  85  The findings of significant differences in the left M l with the 3DMl-baesd method is intriguing, and suggests that the spatial distribution within the M l is different despite the similarity in level of activation. The dorsolateral prefrontal cortex ( D L P F C ) has been shown to be involved in many cognitive-motor tasks, including sequence learning [35], spatial working memory [36], and movement planning and execution based on visual cues [37]. The finding of differential involvement in the D L P F C between E G and 1G movements supports the idea of task-related differences in the functional geographic distribution of activation within this relatively large anatomical area.  3.5.2  Spatial Analysis of f M R I Data  Our results suggest that although the level of activation might have been modulated by task, the changes in the spatial distribution appeared to be much more discriminant. While it may be possible to envisage a sensitive method that merely amplifies inter-subject variability, the proposed method based on invariant spatial descriptors was able to capture task-modulated spatial aspects of activation statistics that are consistent across subjects. Thus, incorporating spatial information could in fact enhance the sensitivity of the activation discrimination analysis. A n important aspect of ROI-based fMRI analysis is the choice of T-threshold. In this thesis, we have decided to use a T-threshold of zero for two reasons. First, there is little consensus on the interpretation of negative T-values [14].  Secondly,  false positives will not likely reside in the same location within an ROI for multiple subjects. Thus, inclusion of these voxels will only decrease discriminability. However, if a high T-threshold is set while some active voxels are buried by noise, these active voxels, which would be consistent across subjects, will be excluded from the  86  analysis. Therefore, using a T-threshold seems to be the most appropriate, and would provide a conservative estimate of the activation differences, while retaining most of the original information of the ROfs. A drawback of the proposed spatial features is that we cannot easily state whether a given region is "more active" or "less active", just that the activation pattern is similar or different. This is analogous to a population-coding perspective at the neuronal level, where stating that a given region is more or less active is very difficult when only the spatial distribution of activation changes. In most cases where there is a significant difference in the spatial distribution of activation statistics, however, the changes in the thresholded mean activation statistics may give some indication of which condition results in "more" activation, even if the difference in the thresholded mean activation statistics does not reach statistical significance due to, e.g. inappropriate selection of threshold.  3.6  Conclusions  In this chapter, we presented a novel approach for fMRI data analysis that incorporates spatial information in contrast to traditional magnitude-based methods. The proposed technique uses 3DMIs to characterize the spatial distribution of brain activations within an ROI. Applying the proposed method to real fMRI data demonstrated that certain spatial aspects of activation statistics appear to be significantly modulated by task performance and being well conserved across subjects. Naively incorporating the spatial information would normally worsen inter-subject variability, but the invariance properties of the proposed spatial descriptors account for variability in subject's orientation in the scanner, while conserving the task-related 'textural' changes in brain activations. Thus as suggested by our results, incorporat87  ing spatial information can greatly enhance the sensitivity in detecting changes in activation statistics, as compared to using mean activation statistics and percentage of activated voxels as features. A direct extension of the proposed method is to look at the spatial distribution changes in time to infer functionally connectivity, a path currently being explored.  88  Bibliography [1] J. L. Lancaster, L. H. Rainey, J. L. Summerlin, C. S. Freitas, P. T. Fox, A. C. Evans, A. W. Toga, and J. C. Mazziotta. Automated labeling of the human brain: A preliminary report on the development and evaluation of a forwardtransform method. Human brain mapping, 5(4):238-242, 1997. [2] A. Nieto-Castanon, S. S. Ghosh, J. A. Tourville, and F. H. Guenther. Region of interest based analysis of functional imaging data. Neurolmage, 19(4): 13031316, 2003. [3] R. T. Constable, P. Skudlarski, E. Mencl, K . R. Pugh, R..K. Fulbright, C. Lacadie, S. E. Shaywitz, and B. A. Shaywitz. Quantifying and comparing regionof-interest activation patterns in functional brain mr imaging: methodology considerations. Magnetic resonance imaging, 16(3):289-300, Apr 1998.  [4] G. Allen and E. Courchesne. Differential effects of developmental cerebellar abnormality on cognitive and motor functions in the cerebellum: an fmri study of autism. The American Journal of Psychiatry, 160(2):262-273, Feb 2003. [5] H. J. Aizenstein, K . A. Clark, M . A. Butters, J. Cochran, V. A. Stenger, C. C. Meltzer, C. F. Reynolds, and C. S. Carter. The bold hemodynamic response in healthy aging. Journal of cognitive neuroscience, 16(5):786-793, 2004. 89  [6] P. Bogorodzki, J. Rogowska, and D. A. Yurgelun-Todd. Structural group classification technique based on regional fmri bold responses. IEEE  Transactions  on Medical Imaging, 24(3):389-398, Mar 2005. [7] B. Ng, R. Abugharbieh, X . Huang, and M . J. McKeown. Characterizing fmri activations within regions of interest (rois) using 3d moment invariants. IEEE Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, NY, June 17-18, 2006, pp. 63 (8 pages). [8] P. Dayan and L. F. Abbott. Theoretical neuroscience: computational and mathematical modeling of neural systems. MIT Press, 2001.  [9] J. V. Haxby, M . I. Gobbini, M . L. Furey, A. Ishai, J. L. Schouten, and P. Pietrini. Distributed and overlapping representations of faces and objects in ventral temporal cortex. Science, 293(5539):2425-2430, 2001. [10] J. F. Mangin, F. Poupon, E. Duchesnay, D. Riviere, A. Cachia, D. L. Collins, A. C. Evans, and J. Regis. Brain morphometry using 3d moment invariants. Medical image analysis, 8(3):187—196, Sep 2004.  [11] K . J. Friston, A. P. Holmes, K. J. Worsley, J. B. Poline, C. D. Frith, and R. S. J. Frackowiak. Statistical parametric maps in functional imaging: a general linear approach. Human Brain Mapping, 2(3): 189-210, 1995. [12] H. Damasio. Human Brain Anatomy in Computerized Images. Oxford Univer-  sity Press, 1995. [13] M . M . Lewis, C. G. Slagle, A. B. Smith, Y . Truong, P. Bai, M . J. McKeown, R. B. Mailman, A. Belger, and X . Huang. Task specific influences of parkinson's  90  disease on the striato-thalamo-cortical and cerebello-thalamo-cortical motor circuitries. Neuroscience, 147(l):224-235, Jun 15 2007. [14] N . Harel, S. P. Lee, T. Nagaoka, D. S. Kim, and S. G. Kim. Origin of negative blood oxygenation level-dependent fmri signals. Journal of cerebral blood flow and metabolism : official journal of the International Society of Cerebral Blood  Flow and Metabolism, 22(8):908-917, 2002. [15] F. A. Sadjadi and E. L. Hall. Three-dimensional moment invariants.  IEEE  Transactions on Pattern Analysis and Machine Intelligence, 2:127-136, 1980.  [16] C. H. Lo and H. S. Don. 3-d moment forms: Their construction and application to object identification and positioning. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(10):1053-1064, 1989.  [17] T. H. Reiss. Features invariant to linear transformations in 2d and 3d. pages 493-496. 11th IAPR International Conference onPattern Recognition, The Hague, Netherlands, 1992. [18] D. Cyganski and J. Orr. Object recognition and orientation determination by tensor methods. Advances in Computer Vision and Image Processing, 3:101-  144, 1988. [19] D. Cyganski and J. Orr. Applications of tensor theory to object recognition and orientation determination. IEEE  Transactions on Pattern Analysis and  Machine Intelligence, 7(6):662-673, 1985.  [20] B. P. Flannery, W. H. Press, S. A. Teukolsky, and W. T. Vetterling. Numerical recipes in c. Press Syndicate of the University of Cambridge, New York, 1992.  91  [21] P. J . Huber. Robust Statistic. Wiley, New York, 1981. [22] R. O. Duda, P. E . Hart, and D. G . Stork.  Pattern  Classification.  Wiley-  Interscience, 2000. [23] J . Ludbrook and H . Dudley. Why permutation tests are superior to t and f tests in biomedical research. The American Statistician, 52(2): 127-132, 1998. [24] S. Vetsa, M . Styner, S. Pizer, J . Lieberman, and G . Gerig. Caudate shape discrimination in schizophrenia using template-free non-parametric tests, pages 661-669. Medical Image Computing and Computer-Assisted InterventionM I C C A I 2003 6th International Conference, Montreal, Canada, Nov 15-18, 2003. [25] R. Benecke, J . C . Rothwell, J . P. Dick, B . L . Day, and C . D. Marsden. Performance of simultaneous movements in patients with parkinson's disease. Brain : a journal of neurology, 109(4):739-757, Aug 1986. [26] R. Benecke, J . C . Rothwell, J . P. Dick, B . L . Day, and C . D. Marsden. Disturbance of sequential movements in patients with parkinson's disease. Brain : a journal of neurology, 110(2):361-379, Apr 1987. [27] K . E . Martin, J . G . Phillips, R. Iansek, and J . L . Bradshaw. Inaccuracy and instability of sequential movements in parkinson's disease. Experimental brain research, 102(1):131-140, 1994. [28] F . Debaere, N . Wenderoth, S. Sunaert, P. Van Hecke, and S. P. Swinnen. Internal vs external generation of movements: differential neural pathways involved in bimanual coordination performed in the presence or absence of augmented visual feedback. Neurolmage, 19(3):764-776, Jul 2003. 92  [29] A . Solodkin, P. Hlustik, D. C . Noll, and S. L . Small. Lateralization of motor circuits and handedness during finger movements. European journal of neurology : the official journal of the European Federation of Neurological Societies, 8(5):425-434, Sep 2001. [30] R. Chen, L . G . Cohen, and M . Hallett. Role of the ipsilateral motor cortex in voluntary movement. The Canadian journal of neurological sciences, 24(4):284291, Nov 1997. [31] R. Chen, C . Gerloff, M . Hallett, and L . G . Cohen. Involvement of the ipsilateral motor cortex in finger movements of different complexities. Annals of Neurology, 41(2):247-254, Feb 1997. [32] J . M . Orgogozo and B . Larsen. Activation of the supplementary motor area during voluntary movement in man suggests it works as a supramotor area. Science, 206(4420):847-850, Nov 16, 1979. [33] P. E . Roland, B . Larsen, N . A . Lassen, and E . Skinhoj. Supplementary motor area and other cortical areas in organization of voluntary movements in man. Journal of neurophysiology, 43(1) :118—136, Jan 1980. [34] G . Goldberg. Supplementary motor area structure and function: review and hypotheses. Behavioral and brain sciences, 8(4):567-615, 1985. [35] J . B. Pochon, R. Levy, J . B. Poline, S. Crozier, S. Lehericy, B. Pillon, B. Deweer, D. Le Bihan, and B. Dubois. The role of dorsolateral prefrontal cortex in the preparation of forthcoming actions: an fmri study. Cerebral cortex (New York, N.Y. : 1991), ll(3):260-266, Mar 2001.  93  [36] G. McCarthy, A. Puce, R. T. Constable, J. H. Krystal, J. C. Gore, and P. Goldman-Rakic. Activation of human prefrontal cortex during spatial and nonspatial working memory tasks measured by functional mri. Cerebral cortex (New York, N.Y. : 1991),  6(4):600-611, Jul-Aug 1996.  [37] S. Funahashi, C. J. Bruce, and P. S. Goldman-Rakic. Mnemonic coding of visual space in the monkey's dorsolateral prefrontal cortex. Journal of neurophysiology, 61(2):331-349, Feb 1989.  94  Chapter 4  Joint Spatial Denoising and Active ROI Delineation 4.1  Introduction  One main application of functional magnetic resonance imaging (fMRI) is in making statistical inferences about the brain mechanisms responsible for different behavioral tasks. Conventionally, this entails generating activation statistics maps with each voxel assigned a statistical value related to its probability of being activated. To make group inferences under this approach, spatial warping of each subject's brain to a common exemplar shape is often performed [1]. However, spatial normalization, typically followed by spatial smoothing, may pool responses from functionally dissimilar regions [2], An alternative approach is to manually draw regions of interest (ROIs) individually for each subject, and examine the statistical properties A version of this chapter has been accepted. B. Ng, R. Abugharbieh, S. J. Palmer, and M. J. McKeown, "Joint Spatial Denoising and Active Region of Interest Delineation in Functional Magnetic Resonance Imaging," 29 Annual International Confernce of the IEEE Engineering in Medicine and Biology Society, Lyon, France, Aug. 23-26, 2007 (accepted, 4 pages). 3  th  95  of regional activation across subjects. Though this ROI-based approach deals with inter-subject anatomical variability, it suffers from two main drawbacks. First, delineating the boundaries of the ROIs is typically guided by anatomical structures [3], which ignores functional properties. Secondly, any outlier voxels can skew the statistical analysis. To deal with Type I errors (i.e. voxels falsely detected as activated) and to increase signal to noise ratio (SNR), a fixed Gaussian filter is conventionally used to spatially smooth the activation statistics maps. This procedure assumes that the M R scanner has a Gaussian-like point spread function (PSF) that smears the fMRI signal over a small neighborhood. Thus based on matched filter theorem [4], the S N R will be increased. However, smoothing the activation maps may affect voxels that do not contain Type I errors, and degrade the spatial resolution. To enhance localization, researchers have investigated spatial wavelet transforms as an alternative [5]. This approach involves representing clusters of voxels with a few wavelet coefficients and de-weighting or removing those coefficients that fall under a certain threshold (i.e. soft or hard thresholding). The inverse wavelet transform is then applied to the surviving coefficients to estimate a denoised activation map. Under this framework, better localization is obtained since not all voxels are altered as in the case of Gaussian smoothing. In this chapter, we propose a method for simultaneously denoising activation statistics within an R O I and refining R O I delineation. Assuming that activation statistics within a small neighborhood are spatially correlated, voxels within each neighborhood would thus exhibit similar level of influence on the overall ROI's response. To impose this condition, isolated voxels at region boundaries with disproportional influence on a derived R O I feature are removed, thereby refining the R O I delineation. Other overly influential voxels are de-weighted based on their influence relative to their neighbors to reduce the  96  biasing effects of these voxels in later analysis. Thus, only voxels detected as exceedingly influential are affected, while most of the original information is retained. The validity of the proposed approach is suggested by the fact that using one feature for denoising (e.g. spatial variance) resulted in greater effect size in another feature (e.g. average activation statistics magnitude).  4.2 4.2.1  Data Acquisition Study Subjects and Experimental Conditions  After obtaining informed consent, fMRI data were collected from 8 healthy subjects. Each subject was required to perform a right-handed motor task that involved squeezing a bulb with sufficient pressure such that an 'inflatable ring', shown as a black horizontal bar on a screen, was kept within an undulating pathway. The pathway is straight during rest periods and becomes sinusoidal at time of stimulus. Each run lasted 260 s, alternating between rest and stimulus of 20 s duration. At time of stimulus, the subject was required to squeeze the bulb at three different frequencies: 0.25 Hz (slow), 0.5 Hz (medium) or 0.75 Hz (fast). The data were collected as part of a larger experiment exploring the rate of change of force production in older subjects and subjects with Parkinson's disease.  4.2.2  Functional M R Imaging  Functional M R I was performed on a Philips Gyroscan Intera 3.0 T scanner (Philips, Best, Netherlands) equipped with a head-coil. We collected echo-planar (EPI) T2*weighted images with B O L D contrast. Scanning parameters were: repetition time 1985 ms, echo time 3.7 ms, flip angle 90°, field of view (FOV) 216x143x240 mm, in  97  plane resolution 128x128 pixels, pixel size 1.9x1.9 mm. Each functional run lasted 4 minutes where 36 axial slices of 3 mm thickness were collected in each volume, with a gap thickness of 1 mm. We selected slices to cover the dorsal surface of the brain and included the cerebellum ventrally. A high resolution 3D Tl-weighted image consisting of 170 axial slices was acquired of the whole brain to facilitate anatomical localization of activation for each subject.  4.2.3  Preprocessing and Activation Statistics Generation  The fMRI data was pre-processed for each subject, using Brain Voyager's (Brain Innovation B . V . ) trilinear interpolation for 3D motion correction and sine interpolation for slice time correction. Further motion correction was performed using motion corrected independent component analysis ( M C I C A ) [6]. No temporal or spatial smoothing was performed on the data. Using Amira (Mercury Computer Systems, San Diego, U S A ) , custom scripts to co-register the anatomical and functional images were generated. Ten ROIs were manually drawn on each structural scan based upon anatomical landmarks and guided by a neurological atlas [3]. ROIs chosen included the thalamus, cerebellum, primary motor cortex ( M l ) , supplementary motor area (SMA), and prefrontal cortex (PFC). The labels on the segmented anatomical scans were resliced at the fMRI resolution. A hybrid Independent Component Analysis (ICA) General Linear Model ( G L M ) scheme [7] was used to contrast each of the three frequency blocks with the periods when subjects were producing static force in generating each subject's activation statistics map. The underlying idea in this scheme is that the G L M regressors (reference and confounds) are derived using ICA, and leave-one-out validation  98  is used to select the best model.  4.3  Methods  The goal of our proposed method is twofold; denoising a subject's ROI activation map (hereafter referred to as T-maps) to reduce the biasing effects of voxels deemed falsely active as well as refining the ROI delineation. Assuming that the ROI T maps have been calculated for two groups (e.g. left S M A for subjects performing two different motor tasks), the proposed method proceeds as follows (Fig. 4.1):  ROI T-map  dFeature  Outlier V o x e l s Detection  Outlier V o x e l s De-weighting  Denoised T-map  Figure 4.1: Framework of proposed spatial denoising method  For each subject's ROI T-map, we first calculate the influence of each voxel on the value of a chosen ROI feature. Average magnitude, percentage of activation, and spatial variance are examined in this study. We then use a robust outlier detection method to determine which voxel(s) appear to be overly 'influential.' Each outlier voxel is then de-weighted or removed by examining its relative influence as compared to its neighbors.  4.3.1  ROI Feature Extraction  Three features are examined in this thesis, namely average magnitude, percentage of activation, and spatial variance. We chose these features, because average magnitude and percentage of activation are widely used in ROI-based fMRI analyses  99  [8], and spatial variance has recently been shown to provide enhanced sensitivity in characterizing the activation distribution within an ROI [9]. Average magnitude is calculated by averaging voxel T-values greater than or equal to zero over the ROI. Percentage of activation is defined as the number of voxels with T-values above 1.96 divided by the total number of voxels within the ROI. Spatial variance, J\, which is based on 3D moments, is calculated as follows:  Jl = m00 + V020 + 77002,  (4.1)  (4.2) ^000 oo  /  poo poo  I  I  (x - x) {y -y) (z p  q  - z) p(x,y,z)dxdydz, r  (4.3)  -oo J —oo J—oo  where n = p+q+r is the order of the centralized 3D moment, p , pqr  (x,y,z) are the  coordinates of a voxel, p(x,y,z) is the T-value of a voxel located at (x,y,z), and x, y, and z are the centroid coordinates of p(x,y,z).  We note that J\ is invariant to  translation, rotation, and scaling, which accounts for differences in brain size across subjects and ensures that the calculated feature values does not depend on subject's orientation in the scanner during acquisition [9]. The invariance properties, thus, allow the 3\ values to be easily compared across subjects.  4.3.2  Voxel Influence Map Generation  To generate the voxel influence map of an ROI for a given subject, we first calculate the feature value, F u, using all voxels in that subject's T-map. We then remove a  a voxel from the original T-map, and recalculate the ROf feature, Fi, with that voxel removed. The 'influence' of the removed voxel is then defined as the absolute difference between F u and Fi, which we denote by dFj. The voxel influence map a  100  of a subject's ROI is then generated by performing this procedure for every voxel within the ROI.  4.3.3  Outlier Voxel Detection  To detect outliers, i.e. overly influential voxels within a subject's ROI, a robust criterion based on the median and the median of absolute deviations ( M A D ) of dFi [10] is used:  lout = find{\ dFi - dF  med  |) > N  std  • a,  (4.4)  where I^t is the set of voxels detected as outliers relative to the other voxels in the ROI, dF  med  dF  med  is the median of dFi, and N t s  d  is the number of standard deviations from  beyond which the value is considered an outlier. Assuming dF is normallyt  distributed, we set N  std  = 3 to retain 99.7% of the possible dFi values as non-outliers.  a is a robust estimate of the standard deviation of dFi, defined as:  a as A • median(\ dFi — dF  med  |),  (4-5)  where A converts M A D to standard deviation and is set to 1.483 assuming dFi is normally-distributed [10]. It is worth noting that our defined criterion (4.4) is robust to outliers, since the median is not as sensitive to extreme values as the mean.  If dF  mean  is in-  stead used in (4.5), extreme outliers can cause a to become arbitrarily large, thus intermediate outliers will likely be undetected.  4.3.4  Outliner Voxel De-weighting  To deweight the T-value (Tk) of each outlier voxel k e lout within a subject's ROI, we first separate its 26-connected neighbors (including itself) into two groups based  101  on (4): outliers  (Vout)  and non-outliers  The number of neighbors may  {V -out)non  be less than 26, especially if voxel k lies at the ROI boundary. Also, errors in ROI delineation can skew the analysis, thus outliers located at region boundaries are treated more stringently than those residing inside the ROI. If voxel k lies inside the ROI and has only outlier neighbors, we leave T unk  modified, since a group of voxels in the same neighborhood exhibiting huge influence on an ROI's feature are more likely to correspond to ROI response. In contrast, if some of voxel k's neighbors are non-outliers, we deweight Tk as follows: T  = %non-out  k  j&Vnon-out  W  where  % n-out n0  j  and  =  1  % t ou  /  d  F  >  WjTj + Voout ^ meVout  l/dF Wm  = ^ ^ - l / d F  jt  >  (-)  ,  (4.7)  WmTm  m  4  6  are the percentage of non-outlier and outlier neighbors.  Our choice of w is based on the assumption that activation statistics within a small neighborhood are spatially correlated. Therefore, no single voxel should exhibit exceedingly greater influence than its neighbors. To impose this condition, we scale down  Ti  by 1/dFi, which proportionally reduces the influence of voxel  i.  Hence,  each of the two sums in (4.6) corresponds to a representative T-value estimate of the non-outlier and outlier voxels. Tk is then approximated by summing the two T-value estimates weighted by  % on-out n  and  %  o  u  t ,  since if %  o u t  »  %non-<mt, Tk is  more likely to correspond to ROI response and vice versa. If voxel k lies on the ROI boundary, Tk is set to zero, except if all its neighbors are outliers, where we deweight Tk using (4.6). Zeroing the T-values in effect refines the ROI boundary, thus accounts for slight errors in ROI delineation.  102  4.3.5  Activation Discriminant Analysis  In this study, we are interested in detecting frequency-related activation differences (fast versus slow frequency) within various ROfs. Given either the raw or the denoised T-maps, we first calculate the average magnitude, T , percentage of actiave  vation, %Act)  ^ d spatial variance, J%, and test whether the differences in feature n  values between fast and slow frequencies, which we denote as activation map feature differences ( A M F D ) , are significant using a paired T-test. The significant level is set at a — 0.05, which corresponds to | T u | of 2.36 for 8 subjects. cr  4.3.6  Comparison of Spatial Denoising Methods  With raw ROf T-maps as baseline, we compared our proposed method to the effects of traditional Gaussian spatial smoothing with a kernel size of 2 voxels (6 mm) full width half max ( F W H M ) . In our approach, the feature used for denoising can be different from that used for activation discriminant analysis. We therefore compared the effects of using T  ave  versus J\ as denoising feature. Also, using different features  for denoising and activation discriminant analysis serves as a validity check. Note that %Act is ° t n  a  suitable denoising feature, since a voxel is either greater or less  than 1.96, thus the resulting influence map will have only two possible values.  4.4  Results and Discussion  Table 4.1 summarizes the activation discriminant analysis results. In general, T  ave  and %Act  increased with task-frequency (i.e. positive T-values), and J\ decreased  with increasing frequency (i.e. negative T-values), thus demonstrating a 'focusing effect'. With the raw ROI T-maps, significant A M F D were detected only in the left  103  Table 4.1: T-values of Activation Map Feature Differences ( A M F D ) Comparing Fast versus Slow Frequency Methods  Proposed Method  Gaussian Smooth  Raw  Test Denoise  Features  Tave  %Act  Jl  Tave.  %Act  Jl  LTHA RTHA LCER RCER LMl RMI LSMA RSMA LPFC RPFC  1.27 1.33 1.18 1.35 1.37 1.80 1.97 1.64 2.43* 2.19  2.15 2.28 1.47 2.03 1.62 2.20 2.61* 1.49 1.90 1.99  -1.64 -1.33 -0.50 -1.16 -2.17 -1.89 -2.41* -1.17 -1.54 -19.3  1.33 1.50 1.55 1.55 1.49 2.00 2.15 1.71 2.66* 2.82*  2.00 1.38 1.59 1.98 1.74 1.82 2.49* 1.76 1.70 2.22  -1.98 -1.74 -1.02 -1.03 -2.22 -2.09 -2.44* -1.07 -1.78 -1.76  THA  Tave Tave Jl  1.28 1.65 1.39 1.56 1.62 1.66 2.23 1.54 3.20* 2.58*  = thalamus, C E R = cerbellum, L = left, R = right, T e aV  Jl  %Act  1.41 1.90 1.38 1.61 2.11 2.01 2.47* 1.52 3.46* 2.22  1 ave 1.91 2.73* 1.34 2.03 1.61 2.27 2.73* 1.58 3.24* 2.52*  Jl  I ave  Jl  2.44* 3.45* 1.37 1.99 1.78 2.86* 3.44* 1.54 3.46* 2.31  -1.14 -2.27 -0.15 -1.39 -1.71 -2.24 -1.66 -1.59 -2.74* -1.71  -2.42* -2.56* -0.83 -1.30 -2.41* -2.54* -2.77* -1.75 -2.64* -1.46  = average magnitude, %Act = percentage of  activation, J\ = spatial variance. First row under "Proposed Method" = feature being tested, second row = feature used for denoising. ^significant at | T it cr  P F C (with T ) ave  | = 2.36.  and left S M A (with %Act and J\). Spatially smoothing the T-maps  with a Gaussian kernel only additionally detected the right P F C (with T ).  In  ave  contrast, the proposed method detected significant A M F D in substantially more ROIs, including the bilateral thalamus, M l , P F C and the left S M A . Activation in the thalamus has been shown to scale with the rate of change of force production [11], and increasing movement speed is known to recruit widespread cortical areas including the right M l [12] and the S M A . Thus, the results generated using the proposed method, are quite consistent with previous neuroscience studies. Comparing T  ave  and J\ as denoising feature, though T  ave  led to some increase  in discriminability (i.e. increase in the T-value magnitude, especially when %Act  w  a  s  used as the test feature) by merely modifying 15% of the voxels within each ROI on average (5% removed, 10% deweighted), substantially greater increase in discriminability was gained by using J\ as the denoising feature (first and third column from the right in Table 4.1), which modified about 20% of the voxels within each ROI (10% removed, 10% deweighted). The reason for the greater discriminability increase is that J\ accounts for voxel location, whereas T  ave  104  would modify any voxel  Figure 4.2: Left thalamus denoised with J\. Shaded region = outlier voxels (~20% of the ROI). Most outlier voxels (~80%) lie on the ROI boundary.  within an ROI regardless of its location. In fact, using J\ as denoising feature, we observed that on average, 80% of the outlier voxels were at the ROI boundary (e.g. Fig. 4.2), compared to 70% when T  ave  was used. This illustrates the importance  of accurate ROI delineation in ROI-based fMRI analysis. Also, we note that using one feature for denoising resulted in greater discriminability in another feature, thus validating our proposed method.  4.5  Conclusion  We presented a new method for refining ROI delineation and denoising fMRI activation statistics within an ROI. Assuming voxels inside a small neighborhood are spatially correlated, we deweight only those voxels that exhibit exceedingly greater influence on an ROI's feature.  Applying our method to real data demonstrated  enhanced discriminability in detecting subtle activation differences as compared to  105  Gaussian spatial smoothing. A direct extension would be to apply our method to raw fMRI signals, which could spatially and temporally denoise the data.  106  Bibliography [1] J . L . Lancaster, L . H . Rainey, J . L . Summerlin, C . S. Freitas, P. T . Fox, A . C . Evans, A . W . Toga, and J . C . Mazziotta. Automated labeling of the human brain: A preliminary report on the development and evaluation of a forwardtransform method. Human brain mapping, 5(4):238-242, 1997. [2] A . Nieto-Castanon, S. S. Ghosh, J . A . Tourville, and F . H . Guenther. Region of interest based analysis of functional imaging data. Neurolmage, 19(4):13031316, 2003. [3] J . Talairach and P. Tournoux. Co-planar stereotaxic atlas of the human brain: an approach to medical cerebral imaging. Stuttgart, Germany: Thieme, 1988. [4] K . J . Worsley, S. Marrett, P. Neelin, and A . C . Evans. Searching scale space for activation in pet images. Human brain mapping, 4(l):74-90, 1996. [5] D . Van De Ville, T . Blu, and M . Unser.  Integrated wavelet processing and  spatial statistical testing of fmri data. Neurolmage, 23(4):1472-1485, Dec 2004. [6] R. Liao, J . L . Krolik, and M . J . McKeown. A n information-theoretic criterion for intrasubject alignment of fmri time series: motion corrected independent  107  component analysis. IEEE Transactions on Medical Imaging, 24(l):29-44, Jan 2005. [7] M . J . McKeown. Detection of consistently task-related activations in fmri data with hybrid independent component analysis. Neurolmage, ll(l):24-35, 2000. [8] R. T . Constable, P. Skudlarski, E . Mencl, K . R. Pugh, R. K . Fulbright, C. Lacadie, S. E . Shaywitz, and B. A . Shaywitz. Quantifying and comparing regionof-interest activation patterns in functional brain mr imaging:  methodology  considerations. Magnetic resonance imaging, 16(3):289-300, Apr 1998. [9] B. Ng, R. Abugharbieh, X . Huang, and M . J . McKeown. Characterizing fmri activations within regions of interest (rois) using 3d moment invariants.  IEEE  Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, N Y , June 17-18, 2006, pp. 63 (8 pages). [10] P. J. Huber. Robust Statistic. Wiley, New York, 1981. [11] D. E . Vaillancourt, M . A . Mayka, K . R. Thulborn, and D . M . Corcos. Subthalamic nucleus and internal globus pallidus scale with the rate of change of force production in humans. Neurolmage, 23(1):175—186, Sep 2004. [12] T . Verstynen, J . Diedrichsen, N. Albert, P. Aparicio, and R. B. Ivry. Ipsilateral motor cortex activity during unimanual hand movements relates to task complexity. Journal of neurophysiology, 93(3):1209-1222, 2005.  108  Chapter 5  Temporal Dynamics of Spatial Distributions in fMRI B O L D Signals 5.1  Introduction  The most common application of functional magnetic resonance imaging (fMRI) is in mapping neural region(s) to particular function(s) by examining which brain areas activate when a certain task is performed. Most conventional analysis methods, such as statistical parametric mapping (SPM) [1], analyze each voxel's timecourse independently and assign a statistics value to that voxel based on its probability of being activated. To make group inferences under this approach, spatial warping of each subject's brain to a common exemplar shape is often performed to create a correA version of this chapter has been accepted. B. Ng, R. Abugharbieh, S. J. Palmer, and M . J. McKeown, "Characterizing Task-Related Temporal Dynamics of Spatial Activation Distributions in fMRI B O L D Signals," 10 International Conference on Medical Image Computing and Computer Assisted Intervention, Brisbane, Australia, Oct. 29-Nov. 2, 2007 (accepted, 8 pages). 4  th  109  spondence between voxels across subjects [2], However, spatial n o r m a l i z a t i o n , w h i c h is t y p i c a l l y followed by spatial smoothing, may inappropriately p o o l responses from functionally dissimilar regions [3], thus degrading i m p o r t a n t spatial information. A n alternative approach that involves d r a w i n g regions of interest (ROIs) i n d i v i d u ally for each subject, and e x a m i n i n g the statistical properties of regional a c t i v a t i o n across subjects, has been shown to offer finer l o c a l i z a t i o n and increased sensitivity to task-related effects [3]. T h i s subject-specific R O I - b a s e d approach is thus followed i n this study. T o determine whether a n R O I is activated or not, a simple approach is to calculate the average intensity over an R O I at every time point, and  determine  if the resulting average time course significantly correlates w i t h the stimulus [4]. T h i s approach, however, ignores any spatial information of a c t i v i t y w i t h i n an R O I and assumes that only signal a m p l i t u d e is m o d u l a t e d by task.  However, spatial  information might be an i m p o r t a n t a t t r i b u t e of b r a i n activity. P r e l i m i n a r y evidence s u p p o r t i n g this idea of spatial characterization was first shown by T h i c k b r o o m et a l . [5], where the spatial extent of activation, as opposed to response magnitude,  was  found to be m o d u l a t e d by different levels of force d u r i n g a sustained finger flexion task.  T h e i r results were based on v i s u a l inspection and counting the number of  activated voxels w i t h i n an R O I . Recently, we presented a more elaborate study of the spatial patterns of a c t i v i t y w i t h i n an R O I where quantitative measures of invariant spatial properties were used to discriminate task-related differences i n b r a i n a c t i v i t y [6]. Results demonstrated that, by e x a m i n i n g changes i n different spatial aspects of an activation d i s t r i b u t i o n , sensitivity i n detecting functional changes is enhanced as compared to using intensity means only. P r e v i o u s analyses e x a m i n i n g spatial patterns of activation, i n c l u d i n g that i n  110  [6], were performed on T-maps where the spatial information is collapsed over time, thus only considered the time-averaged spatial patterns of brain activity. In this chapter, we extend our previously proposed spatial characterization approach to the temporal domain to explore whether the spatial distribution of the blood oxygenation level-dependent ( B O L D ) signal itself is modulated in time by task performance. We note an important difference between our current and previous work [6] is that the generated spatial feature time courses can be used to infer R O I activation, as opposed to only comparing 2 groups of time-averaged activation statistical maps. To characterize the spatial changes, three dimensional (3D) moment descriptors were used as features and were calculated at each time point. The magnitudes of these features, however, are normally not comparable across subjects due to inter-subject variability in brain shapes and sizes, but are comparable for the same subject over time. Thus, any detected modulations of the spatial features over time for a given subject may in fact represent meaningful spatial changes in activation. In this study, eight healthy subjects were recruited to perform a bulb-squeezing task at various frequencies. The cerebellum, primary motor cortex ( M l ) , supplementary motor area (SMA), prefrontal cortex ( P F C ) , and anterior cingulate cortex ( A C C ) were chosen as regions of interest. We demonstrate that our method can both detect activation within an R O I , as well as discriminate differences in activation patterns at the various task frequencies. This confirms previous findings of the value in incorporating spatial information into traditional intensity-based fMRI analyses.  Ill  Figure 5.1: Experimental task and stimulus timing, (a) Subjects were required to keep the side of the black ring on the gray path (see text), (b) R = rest, Slow, Med, and Fast = stimulus at 0.25, 0.5, and 0.75 Hz, respectively. Each block is 20 s in duration.  5.2  Data Acquisition and Preprocessing  In this study, after informed consent was obtained, fMRf data were collected from 8 healthy subjects. Each subject was required to perform a right-handed motor task that involved squeezing a bulb with sufficient pressure such that an 'inflatable ring', shown as a black horizontal bar on a screen, was kept within an undulating pathway (Fig. 5.1(a)). The pathway remains straight during rest periods and becomes sinusoidal at time of stimulus. Each run lasted 260 s, consisting of a 20 s rest period at the beginning and end, 6 stimuli of 20 s duration, and 20 s rest periods between the stimuli, as shown in Fig. 5.1(b). At time of stimulus, the subject was required to squeeze the bulb at 0.25, 0.5 or 0.75 Hz, corresponding to 'Slow', 'Med', and 'Fast' in Fig. 5.1(b). The data were collected as part of a larger experiment exploring the rate of change of force production in older subjects and subjects with Parkinson's disease.  112  5.2.1  fMRI Data Acquisition  Functional M R I was performed on a Philips Gyroscan Intera 3.0 T scanner (Philips, Best, Netherlands) equipped with a head-coil. We collected echo-planar (EPI) T 2 * weighted images with B O L D contrast. Scanning parameters were: repetition time 1985 ms, echo time 3.7 ms, flip angle 90°, field of view (FOV) 216x143x240 mm, in plane resolution 128x128 pixels, pixel size 1.9x1.9 mm. Each functional run lasted 4 minutes where 36 axial slices of 3 mm thickness were collected in each volume, with a gap thickness of 1 mm. We selected slices to cover the dorsal surface of the brain and included the cerebellum ventrally. A high resolution 3D Tl-weighted image consisting of 170 axial slices was acquired of the whole brain to facilitate anatomical localization of activation for each subject.  5.2.2  fMRI Preprocessing  The fMRI data was preprocessed for each subject, using Brain Voyager's (Brain Innovation B . V . ) trilinear interpolation for 3D motion correction and sine interpolation for slice time correction. Further motion correction was performed using motion corrected independent component analysis (MCICA) [7]. To handle temporal autocorrelations, a 'coloring' scheme was used [8], where the time series were high-pass filtered at 0.02 Hz (task-frequency being 0.025 Hz) to remove the majority of the low frequency noise, and temporally smoothened with a Gaussian of width 2.8s [8]. The first and last 20 s of the time series were truncated to mitigate transient effects. No spatial smoothing was performed. The Brain Extraction Tool ( B E T ) in MRIcro [9] was used to strip the skull off of the anatomical and first functional image from each run to enable a more accurate alignment of the functional and anatomical scans. Custom scripts to co-register the  113  anatomical and functional images were generated using the Amira software (Mercury Computer Systems, San Diego, U S A ) . Ten specific ROIs were manually drawn on each unwarped structural scan using Amira. The following ROIs were drawn separately in each hemisphere, based upon anatomical landmarks and guided by a neurological atlas [10]: cerebellum, M l (Brodman Area 4), S M A (Brodman Area 6), P F C (Brodman Area 9 and 10), and A C C (Brodman Area 28 and 32). The labels on the segmented anatomical scans were resliced at the fMRI resolution. The raw time courses of the voxels within each ROf were then extracted for analysis as described in the next section.  5.3 Methods The main goal of the proposed method is to demonstrate that temporal dynamics of the spatial distribution in B O L D signals can be used to infer whether an ROf is activated, as well as to discriminate differences in activation patterns at various task frequencies. Details of feature time course extraction, activation detection, and activation pattern discrimination are discussed below.  5.3.1  Feature Time Course Extraction  The spatial feature descriptors used in this thesis are based on centralized 3D moments, defined as:  Hvqr(t)=  /  /  (x - x) {y-y) (z p  q  - z) p(x,y,z,t)dxdydz, r  (5.1)  where n = p+q+r is the order of the moment, (x,y,z) are the coordinates of a voxel, p(x, y, z, t) is the intensity of a voxel located at (x,y,z) inside a given ROI at time t, and x, y, and z are the centroid coordinates of p(x,y, z,t). 114  To untangle the effect  of amplitude changes, p(x, y, z, t) is normalized such that the intensity values of the voxels within the ROI sums up to one at every time point t. This step ensures that the mean ROI intensity does not change with time. Thus, any detected modulations in the spatial feature will be purely due to spatial changes in the B O L D signal. To ease interpretation of the results and since higher order moments are less robust to noise [11], only 2  nd  and 3  r d  order 3D moment descriptors characterizing spatial  variance [12] and skewness, respectively, were used:  (5.2)  Jl{t) = M200(*) + M020(*) + M002(*)>  S(t)  =  P30o{t)  +  £i()30(£) +  (5.3)  P003(t)  To compare with the results obtained using the proposed spatial feature time courses, the traditionally used mean intensity time course, I(t), for each ROI of a given subject is calculated by averaging the intensity over the ROI at every time point.  5.3.2  Activation Detection  To make group inference as to whether a given ROI is activated, each subject's ROI feature time courses (spatial or mean intensity) are first correlated with a box-car that is time-locked to stimulus with a delay of 4 s [13]. We did not convolve the box-car with a haemodynamic response function since spatial changes, as governed by the different onsets of the active voxels, may exhibit a different temporal profile than that of the haemodynamic response.  For each subject, this results in thirty  correlation values, one per combination of feature and ROI (e.g.  J\(t), left M l ) .  Each correlation value is then converted into a T-value (5.4):  (5.4)  115  where r is the correlation value and N is the number of samples used in generating r. The set of T-values of a particular combination of feature and ROI from all subjects is then tested against 1.96 using a T-test to determine the probability (p-value) that the T-values are lower than 1.96 (i.e. the probability that ROI is not activated). The critical p-value was chosen at 0.05.  5.3.3  Activation Pattern Discrimination  To discriminate the differences in activation pattern at the various task frequencies, each subject's ROI feature time courses (spatial or mean intensity) are first segmented according to Fig. 5.2. Except for the first and last segments, each segment consists of a 10 s rest before and after the 20 s stimulus. Segments of the same task frequency are concatenated and correlated with the corresponding segments of the shifted reference signal (see Fig. 5.2). This results in ninety correlation values per subject, one for each combination of frequency, feature, and ROI (e.g. slow,  J\(t),  left M l ) . Each correlation value is then converted into a T-value using (5.4). For each combination of feature and ROI, the set of T-values of a particular frequency from all subjects are tested pair-wise against the other two frequencies (i.e. fast versus slow, fast versus medium, medium versus slow). This is performed using a T-test to determine the probability (p-value) that the sets of T-values from the two task frequencies are the same (i.e. the probability that activation patterns at the two frequencies are the same). The critical p-value was chosen at 0.05.  5.4  Results and Discussion  Table 5.1 summarizes the activation detection results generated by extracting the spatial and intensity features from real fMRI data as described in Section 5.3.1, and 116  !_ 1  J  a  I  0  L .  r  1  r  r  r  !  H f ?  K 20  1  r  1  40  i  _  l  60  /  /\ i;  v.  /;  1  1  1  80  100  120 time  U l  140  (sec)  1  160  - I  180  1  200  _  A'  / l  220  1  240  \ 1  260  Figure 5.2: Feature time course segmentation. The box-car curve corresponds to timing of the stimulus delayed by 4 s. The solid line is a sample feature time course (spatial variance, Ji(t), of the left M l averaged over subjects with its temporal mean removed and divided by its standard deviation). The dotted lines show how the feature time courses are parsed into 6 segments. Slow, Med, and Fast correspond to the task frequencies of 0.25, 0.5, and 0.75 Hz, respectively.  correlating the resulting feature time courses with the reference signal. Using spatial variance, Ji(t), the left cerebellum, left M l , both SMAs, left P F C , and left A C C were detected as active. We expected the left M l to be activated, as typically observed for right-handed motor tasks. It is worth noting that the left M i ' s B O L D signal distribution shown reduced spatial variance (i.e. focuses) during the time of stimulus (see Fig. 5.2). Skewness, S(t), additionally detected activation in the right M l . These results demonstrate that the spatial distribution of B O L D signals is, in fact, modulated by stimulus, which supports our hypothesis that spatial changes in B O L D signals are task-related and can be used to infer activation within an ROI. Using the traditional mean intensity measure, the right cerebellum, both SMAs, left P F C , and both A C C s were detected as active. Comparing to the results generated with the proposed spatial features, some consistencies are shown. In fact, based on the results in Table 5.1, activation within an ROI appears to modulate  117  Table 5.1: p-values of ROI Activation Feature Mt) I(t) S(t)  C E R = cerebellum, Ji(t)  LCER  0.002*  0.077  0.469  RCER  0.220  0.066  0.042*  LMl  0.002*  RMI  0.157  0.051 0.047*  0.162 0.164  LSMA  0.027*  0.034*  0.021*  RSMA  0.044*  0.035*  0.034*  LPFC  0.020*  0.034*  RPFC  0.036* 0.341  0.185  0.058  LACC  0.001* 0.056  0.301 0.112  0.023*  RACC  0.045*  = spatial variance, S(t) = skewness, I(t) = mean intensity, L = left, R = right,  'statically significant at a = 0.05.  both in amplitude and in space. Segmenting the feature time courses according to task frequencies and using the proposed spatial features, significant frequency-related activation differences were detected in the right cerebellum and right M l when comparing fast versus slow frequencies (Table 5.2). These results matched our expectations since the modulation of movement speed is known to involve a complex network of brain areas, including the right cerebellum and right M l [14].  Also, significant activation differences were found  in the right P F C and right A C C using the proposed spatial features. In contrast, no significant activation differences were found using mean intensity. Also, no significant activation differences were detected when comparing fast versus medium frequencies and medium versus slow frequencies for any of the features, thus these results were excluded in Table 5.2. Examining the results in Table 5.2, spatial changes appear to provide greater sensitivity in detecting subtle activation differences as compared to intensity. 118  Table 5.2: p-values of Activation Differences Comparing Fast versus Slow Frequencies Feature  Ji(t)  S(t)  I(t)  LCER  0.1834  0.1931  0.1528  RCER  0.9993 0.2932  0.0096* 0.4017  0.2757  LMl RMI  0.0442*  0.3223  0.1524*  LSMA  0.6872  0.2499  0.1836  RSMA  0.7887  0.9066  0.0579  LPFC  0.5981  0.3059  0.2760  RPFC  0.0398*  0.3668  0.7416  •LACC RACC  0.2551  0.4029  0.2126  0.0311*  0.1028  0.3248  0.4190  C E R = cerebellum, Ji(t) = spatial variance, S(t) = skewness, I(t) = mean intensity, L = left, R = right, 'statically significant at a — 0.05.  5.5  Conclusion  In this chapter, we proposed using 3D moment-based spatial descriptors to characterize the temporal dynamics of spatial activation distribution within an ROI for fMRI analysis. We demonstrated with real fMRf data that certain spatial aspects of activation, as opposed to just amplitude, are modulated by stimulus - a result that appeared consistent across subjects. Furthermore, we showed that our method was able to better discriminate frequency-related differences in activation patterns during motor task performance when compared to using mean intensity only. These results suggest that spatial characterization of B O L D signal can complement traditional intensity-based fMRI analysis. A direct extension of the proposed method would be to examine functional connectivity and phase relations between ROIs using spatial feature time courses, an approach currently being pursued.  119  Bibliography [1] K . J . Friston. Statistical parametric mapping and other analyses of functional imaging data. Brain Mapping, The Methods, pages 363-396, 1996. [2] J . L . Lancaster, L . H . Rainey, J . L . Summerlin, C . S. Freitas, P. T . Fox, A . C . Evans, A . W . Toga, and J . C . Mazziotta. Automated labeling of the human brain: A preliminary report on the development and evaluation of a forwardtransform method. Human brain mapping, 5(4):238-242, 1997. [3] A . Nieto-Castanon, S. S. Ghosh, J . A . Tourville, and F . H . Guenther. Region of interest based analysis of functional imaging data. Neurolmage, 19(4): 13031316, 2003. [4] Y . Liu, J H Gao, M . Liotti, Y . Pu, and P T Fox. Temporal dissociation of parallel processing in the human subcortical outputs. Nature, 400(6742):364-7, 1999. [5] G . W . Thickbroom, B . A . Phillips, I. Morris, M . L . Byrnes, and F . L . Mastaglia. Isometric force-related activity in sensorimotor cortex measured with functional  mri. Experimental Brain Research, 121(l):59-64, Jul 1998. [6] B . Ng, R. Abugharbieh, X . Huang, and M . J . McKeown. Characterizing fmri activations within regions of interest (rois) using 3d moment invariants. I E E E  120  Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, NY, June 17-18, 2006, pp. 63 (8 pages). [7] R. Liao, J. L. Krolik, and M. J. McKeown. An information-theoretic criterion for intrasubject alignment of fmri time series: motion corrected independent component analysis. IEEE Transactions on Medical Imaging, 24(l):29-44, Jan 2005. [8] K. J. Friston, A. P. Holmes, J. B. Poline, P. J. Grasby, S. C. Williams, R. S. Frackowiak, and R. Turner. Analysis of fmri time-series revisited. Neurolmage, 2(l):45-53, Mar 1995. [9] C. Rorden. Stereotaxic display of brain lesions.  Behavioural  Neurology,  12(4):191-200, 2000. [10] J. Talairach and P. Tournoux. Co-planar stereotaxic atlas of the human brain: an approach to medical cerebral imaging. Stuttgart, Germany: Thieme, 1988. [11] B. P. Flannery, W. H. Press, S. A. Teukolsky, and W. T. Vetterling. Numerical recipes in c. Press Syndicate of the University of Cambridge, New York, 1992.  [12] F. A. Sadjadi and E. L. Hall. Three-dimensional moment invariants. IEEE Transactions on Pattern Analysis and Machine Intelligence,  2:127-136, 1980.  [13] M. Svensen, F. Kruggel, and DY von Cramon. Probabilistic modeling of singletrial fMRI data. IEEE Transactions on Medical Imaging, 19(l):25-35, 2000. [14] T. Verstynen, J. Diedrichsen, N. Albert, P. Aparicio, and R. B. Ivry. Ipsilateral motor cortex activity during unimanual hand movements relates to task complexity. Journal of neurophysiology, 93(3):1209-1222, 2005.  121  Chapter 6  Conclusions and Future Work A novel ROI-based fMRI analysis paradigm is proposed for quantitatively characterizing the spatial distribution of brain activation within an ROI. We demonstrated for the very first time that the information encoded by voxel locations can be meaningfully incorporated in ROI-based fMRf analysis by using invariant spatial descriptors. The essence of our method is based on the invariant nature of the proposed descriptors, which accounts for inter-subject variability in brain shape, brain size, and subject's orientation in the scanner. Our proposed analysis paradigm facilitates both spatial and spatio-temporal analysis of fMRI data.  6.1  Spatial Analysis of Activation Statistics Maps  Incorporating spatial information into ROI-based fMRI analysis by using invariant spatial descriptors (3DMIs) demonstrated enhanced discriminability of task-related activation differences over using traditional magnitude-based features [1], [2], [3]. However, several important issues remained to be addressed, namely the choice of activation statistics threshold for isolating activated voxels, the multiple comparisons 122  problem, and further validation of the results. In the current method, all voxels with positive T-statistics are used, and the justification is that false positives (voxels falsely deemed active) will not reside at the same location within an ROf across subjects. Hence, inclusion of these voxels will only reduce discriminability. Yet, truly active voxels corrupted by noise may reside in similar location within an ROI across subjects.  Thus, if a "high" T-threshold  is set, some of these active voxels will be excluded from the analysis. Therefore, using a T-threshold of zero provides a conservative estimate of activation differences, while retaining most of the original information of the ROIs. However, developing methods to first isolate the activated voxels before spatial characterization might further enhance discriminability. In fact, the processing method in [3] partially deals with this issue, where overly influential voxels that reside in isolation within an ROI (thus likely to be non-active) are deweighted. By further exploiting the idea in [3], other non-active voxels can be identified. Regarding the multiple comparisons correction issue, all current results are thresholded at an uncorrected p-value of 0.05. The justification is that only tens of ROIs are examined as opposed to thousands of voxels in voxel-based analysis. Also, analogous to the assumption of spatial correlation between voxels in the G R F multiple comparison correction approach, the selected ROIs are all motor related, thus their response might also be correlated. Thereby, using Bonferroni or F D R correction might be too conservative. Instead, perhaps through examining the functional connectivity between ROIs, correction can be applied to networks of ROIs, which accounts for ROI-level spatial correlation, leading to less conservative thresholds. Concerning extended validation of results, all our results were compared to prior neuroscience findings. However, most prior studies only suggest that the  123  magnitude of the ROI responses are modulated by task, whereas we are suggesting that different tasks are encoded by the spatial distribution of activation. Since the ground truth is not known, we are currently relying on indirect supporting evidence from prior studies. To more elaborately validate our method, realistic synthetic data is needed. One way to generate such data is to segment an ROI that is expected to incur no activation, and synthetically inject activation into that ROI. The drawback is that the noise structure across ROIs could be quite different, thus such approach does not guarantee that our results are valid for all ROIs. Alternatively, we can spatially permute the voxel locations and check if significant activation differences are still detected. This procedure is only appropriate for spatial feature analysis, but nonetheless provides a new validation approach to explore.  6.2  Spatio-temporal Analysis of B O L D fMRI Signal  In [4], we demonstrated that the spatial distribution of BOLD signal is modulated in time by task performance, which has significant neuroscience implications since all prior studies assume only the magnitude of BOLD signal is modulated by task. In functional connectivity analysis, an ROI response is often summarized by averaging voxel intensity time courses over an ROI. A more meaningful way to summarize an ROI might be to use spatial feature time course. In fact, examining the phase relations between spatial feature and mean intensity time courses might elucidate other functional properties of the brain. Also, examining the phase relations between ROIs might help infer causal relations. We note that in [4], the issue of removing temporal autocorrelations was not fully addressed, since only the coloring scheme was tested. Comparing results generated by other preprocessing approach like whitening, will be beneficial. Also, 124  we suspect that the temporal autocorrelation problem might be less pronounced in the spatial domain. Thus, comparing the amount of residual autocorrelations in the spatial feature and mean intensity would be interesting. The range of applications of our proposed analysis paradigm is very significant. In fact, with adequate validation, we foresee our approach being prospectively applied in clinical settings for diagnosing neurological diseases and evaluating drug therapies. Also, our paradigm can be easily transferred to examine other neurological diseases due to its generality.  125  B i b l i o g r a p h y  [1] B . Ng, R. Abugharbieh, X . Huang, and M . J . McKeown. Characterizing fmri activations within regions of interest (rois) using 3d moment invariants. f E E E Computer Society Workshop on Mathematical Models in Biomedical Image Analysis, New York, N Y , June 17-18, 2006, pp. 63 (8 pages). [2] B . Ng, R. Abugharbieh, X . Huang, and M . J . McKeown. Characterization of fmri activation maps using invariant 3d moment descriptors. IEEE  Transactions  on Medical Imaging, invited paper, (under revision). [3] B . Ng, R. Abugarbieh, S. J . Palmer, and M . J. McKeown. Joint spatial denoising and active region of interest delineation in functional magnetic resonance imaging.  29  th  Annual International Confernce of the I E E E Engineering in Medicine  and Biology Society, Lyon, France, Aug. 23-26, 2007 (accepted, 4 pages).  [4] B . Ng, R. Abugharbieh, S. J . Palmer, and M . J . McKeown.  Characterizing  task-related temporal dynamics of spatial activation distributions in fmri bold signals. 10  th  International Conference on Medical Image Computing and Com-  puter Assisted Intervention, Brisbane, Australia, Oct. 29-Nov. 2, 2007 (accepted, 8 pages).  126  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0100697/manifest

Comment

Related Items