You may notice some images loading slow across the Open Collections website. Thank you for your patience as we rebuild the cache to make images load faster.

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A robust strategy for cleaning motion artifacts in resting state fMRI Yu, Tianze 2019

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

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
24-ubc_2019_september_yu_tianze.pdf [ 14.12MB ]
Metadata
JSON: 24-1.0379472.json
JSON-LD: 24-1.0379472-ld.json
RDF/XML (Pretty): 24-1.0379472-rdf.xml
RDF/JSON: 24-1.0379472-rdf.json
Turtle: 24-1.0379472-turtle.txt
N-Triples: 24-1.0379472-rdf-ntriples.txt
Original Record: 24-1.0379472-source.json
Full Text
24-1.0379472-fulltext.txt
Citation
24-1.0379472.ris

Full Text

A ROBUST STRATEGY FOR CLEANING MOTION ARTIFACTS  IN RESTING STATE FMRI by  Tianze Yu    A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF  MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Electrical and Computer Engineering)   THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver)  June 2019  © Tianze Yu, 2019 ii    The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled:  A ROBUST STRATEGY FOR CLEANING MOTION ARTIFACTS IN RESTING STATE FMRI  submitted by Tianze Yu in partial fulfillment of the requirements for the degree of Master of Applied Science in Electrical and Computer Engineering  Examining Committee: Z. Jane Wang Supervisor  Martin J. McKeown Co-supervisor Rabab K. Ward Supervisory Committee Member  Additional Examiner   Additional Supervisory Committee Members:  Supervisory Committee Member  Supervisory Committee Member iii  Abstract  There is an increasing recognition that different preprocessing approaches for functional magnetic resonance imaging (fMRI) data may have a profound impact on downstream analyses. A critical element of standard preprocessing of fMRI data is motion correction, as head motion during fMRI scanning can induce changes in blood oxygenation level dependent (BOLD) signals that may confound estimations of brain activity, particularly connectivity estimates between brain regions. In this thesis, we propose an approach which explicitly decouples the changes due to brain activation and changes due to motion. Independent Vector Analysis (IVA) is used to determine the optimal combination of basis volumes to match reference volumes. Such an approach, which we call Motion correction with IVA (McIVA), is amenable to wide-scale parallelization. The mutual information between the first volume and all subsequent volumes declined more slowly after preprocessing with McIVA compared to the raw data and other motion correction schemes. Remarkably, since the final motion-correction volume in McIVA is based on a combination of images, we show that McIVA's error can be actually less than interpolation error. What is more, McIVA resulted in the lowest connectivity across a range of spatial distances between regions of interest (ROI) pairs. When a volume is severely corrupted, the IVA fails to converge, providing a principled way to determine if a given time point is too corrupted for recovery. Finally, we assessed the effects of McIVA on two popular denoising methods, aCompCor and ICA-AROMA on resting-state fMRI derived. McIVA resulted in reducing inflated connectivity estimates, while still retaining an adequate degree of freedom. Though the proposed method is based on resting-state fMRI data, it could be applied to task-related fMRI data as well. We conclude that the proposed iv  approach is superior for removing motion-related artifacts and reducing biases in functional connectivity estimates induced by head movement. v  Lay Summary  When people undergo the functional magnetic resonance imaging (fMRI), it is hard for them to stay dormant and relaxed throughout the entire imaging process which usually lasts for several minutes to tens of minutes. Any little movement of the head, for example, the patient moves his head for 0.1 mm or rotates 0.1 degrees, could have a significant effect on the brain connectivity analysis of fMRI signals. Therefore, an essential pre-processing step of fMRI is motion correction and denoising. In this thesis, a new independent vector analysis (IVA)-based motion correction approach for fMRI is proposed.  vi  Preface  The research was jointly initiated by Dr. Z. Jane Wang, Dr. Martin J. McKeown and the thesis author, and the majority work of the research was conducted by the author of this thesis, with valuable suggestions from Dr. Z Jane Wang and Martin J. McKeown. Dr. Saurabh helped a lot on using SPM for fMRI preprocessing in Chapter 3 and also gave some suggestions on applying the mask on evaluating the ROI-based correlation in Chapter 4. Dr. Aiping Liu helped a lot in analyzing the results in Chapter 4 and Chapter 5.  vii  Table of Contents  Abstract ........................................................................................................................................ iii Lay Summary ................................................................................................................................. v Preface ........................................................................................................................................... vi Table of Contents ......................................................................................................................... vii List of Tables .................................................................................................................................. x List of Figures ............................................................................................................................... xi List of Abbreviations ................................................................................................................. xiii Acknowledgements ...................................................................................................................... xv Chapter 1: Introduction ................................................................................................................ 1 1.1 Overview of functional MRI ........................................................................................... 1 1.2 Properties of fMRI data ................................................................................................... 4 1.3 Interference of fMRI ........................................................................................................ 6 1.3.1 Head movement during scanning ............................................................................ 6 1.3.2 Noise interference of devices and physiological activities ...................................... 8 1.3.3 Susceptibility and N/2 artifacts ................................................................................ 8 1.4 Conventional pre-processing pipeline for fMRI data .................................................... 10 1.4.1 Realignment ........................................................................................................... 11 1.4.2 Co-registration ....................................................................................................... 12 1.4.3 Normalization ........................................................................................................ 12 1.4.4 Smoothing .............................................................................................................. 14 1.4.5 Segmentation ......................................................................................................... 15 viii  1.4.6 Denoising ............................................................................................................... 16 1.4.6.1 Realignment regression ..................................................................................... 17 1.4.6.2 Censoring ........................................................................................................... 17 1.4.6.3 aCompCor .......................................................................................................... 18 1.4.6.4 ICA-AROMA .................................................................................................... 19 1.5 Effect of head motion .................................................................................................... 19 1.6 Motivation, objectives and thesis organization ............................................................. 21 Chapter 2: Motion correction with IVA .................................................................................... 25 2.1 Motivation ..................................................................................................................... 25 2.2 Independent Vector Analysis ......................................................................................... 26 2.3 IVA-L algorithm ............................................................................................................ 27 2.4 IVA-G algorithm ........................................................................................................... 29 2.5 Model structure .............................................................................................................. 30 Chapter 3: Evaluation of McIVA ............................................................................................... 34 3.1 Dataset and Preprocessing ............................................................................................. 34 3.1.1 Participants ............................................................................................................ 34 3.1.2 Data acquisition ..................................................................................................... 34 3.1.3 Data preprocessing ................................................................................................ 35 3.2 Measure methods ........................................................................................................... 37 3.2.1 Absolute frame displacement ................................................................................ 37 3.2.2 Mutual information ................................................................................................ 37 3.2.3 Assessment of interpolation effects ....................................................................... 38 3.2.4 Outlier detection .................................................................................................... 39 ix  3.2.5 Functional connectivity correlations ..................................................................... 39 Chapter 4: Results ....................................................................................................................... 42 4.1 Mutual information between first and subsequent volumes .......................................... 42 4.2 Mean mutual information between two volumes .......................................................... 43 4.3 Mean displacement ........................................................................................................ 44 4.4 Interpolation Error ......................................................................................................... 46 4.5 Evaluation of processing mechanism of outliers ........................................................... 47 4.6 Functional connectivity ................................................................................................. 50 4.7 Comparison with common denoising methods ............................................................. 52 4.7.1 Effect of preprocessing .......................................................................................... 52 4.7.2 Comparison with denoising approaches ................................................................ 55 4.7.3 ROI-based functional correlation .......................................................................... 58 Chapter 5: Discussion .................................................................................................................. 60 5.1 Difference between McIVA and MCICA ..................................................................... 61 5.2 Concern about white matter ........................................................................................... 73 5.3 Number of reference volumes ....................................................................................... 76 5.4 Drawbacks ..................................................................................................................... 78 Chapter 6: Conclusion and Future work .................................................................................. 79 6.1 Conclusion and Contribution ......................................................................................... 79 6.2 Future Considerations .................................................................................................... 81 Bibliography ................................................................................................................................. 82  x  List of Tables  Table 3-1. Characteristics of the participants included in the dataset. .......................................... 34 Table 5-1. Voxel-wised functional connectivity based on different motion conditions. .............. 69  xi  List of Figures  Figure 1-1. Workflow of fMRI acquisition process. ....................................................................... 2 Figure 1-2. Data structure of raw fMRI data. .................................................................................. 3 Figure 1-3. Example of structural data and functional data. ........................................................... 3 Figure 1-4. Displacement caused by head motion during fMRI scanning. ..................................... 7 Figure 1-5. Estimation of head motion. ........................................................................................... 8 Figure 1-6. N/2 ghost artifact in fMRI scan. ................................................................................... 9 Figure 1-7. Conventional pre-processing pipeline for fMRI data. ................................................ 10 Figure 1-8. Spatial normalization of fMRI volume. ...................................................................... 13 Figure 1-9. Spatial smoothing with Gaussian kernel. .................................................................... 15 Figure 1-10. Brain segmentation. .................................................................................................. 16 Figure 1-11. Illustrative examples of correlation changes caused by head movement. ................ 20 Figure 2-1. The main structure of the proposed McIVA. .............................................................. 31 Figure 3-1. Preprocessing pipeline for dataset in Table 3-1. ......................................................... 36 Figure 3-2. The timepoint extracted from real fMRI data and artificially corrupted volume. ...... 39 Figure 4-1. Mutual information between the first and subsequent volumes. ................................ 42 Figure 4-2. Mutual information of the original and the motion corrected fMRI data. .................. 44 Figure 4-3. The average displacement of the original and the motion corrected fMRI data. ........ 45 Figure 4-4. Mean absolute error of two approaches. ..................................................................... 46 Figure 4-5. The voxel intensity changes of two approaches. ........................................................ 47 Figure 4-6. Motion estimation of the original and motion corrected fMRI data with SPM. ......... 48 Figure 4-7. Motion estimation of the fMRI data motion corrected with McIVA. ........................ 49 xii  Figure 4-8. Cross-correlation coefficient in McIVA. .................................................................... 49 Figure 4-9. ROI based function correlation under specific motion correction methods. .............. 51 Figure 4-10. The effect of the preprocessing pipeline. .................................................................. 53 Figure 4-11. Functional connectivity distribution of aCompCor, McIVA and ICA-AROMA. .... 57 Figure 4-12. Functional connectivity between the seed ROI and other ROIs. .............................. 59 Figure 5-1. Structure of the IVA Core. .......................................................................................... 63 Figure 5-2. Motion estimation of the fMRI data motion corrected with MCICA and McIVA. .... 64 Figure 5-3. Relationship of H(y) H(x) and the estimated motion in real fMRI data. .................... 65 Figure 5-4. Convergence properties of the cost function of MCICA. ........................................... 66 Figure 5-5. Convergence properties of the cost function of McIVA. ............................................ 67 Figure 5-6. Voxel-wised functional connectivity distribution on subject level. ........................... 68 Figure 5-7. Distribution of functional connectivity of fMRI data on group level. ........................ 70 Figure 5-8. Distribution of subject-level QC measures (a). .......................................................... 71 Figure 5-9. Distribution of subject-level QC measures (b). .......................................................... 72 Figure 5-10. Motion estimation of the original fMRI data of the subject 4. ................................. 73 Figure 5-11. Motion denoising using McIVA combined with nuisance regression. ..................... 74 Figure 5-12. Group-wise ROI-to-ROI correlation of Figure 5-11. ................................................ 76 Figure 5-13. ROI based correlation under different number of reference volumes of McIVA. ... 77  xiii  List of Abbreviations  BSS – Blind source separation CSF – Cerebral Spinal Fluid DOF – Degree of freedom FC – Functional connectivity FD – Frame-wise displacement fMRI – functional Magnetic Resonance Imaging GLM – General linear model GM – Gray Matter HC – Health Control ICA – Independent component analysis IVA – Independent vector analysis MAE – Mean absolute error McIVA – Motion correction with Independent Vector Analysis MI – Mutual information MNI – Montreal Neurological Institute MRI – Magnetic Resonance Imaging PC – Principle component PD – Parkinson’s Disease ROI – Region of Interest RP – Realignment parameter SCV – Source component vector xiv  WM – White Matter xv  Acknowledgements  I want to express my great appreciation to my supervisors, Dr. Z. Jane Wang and Dr. Martin J. McKeown for their help, support, encouragement, and guidance for the research throughout my master study.  I want to thank my friends and lab mates for their help and encouragement.   I own my deepest gratitude to my parents and my wife for their love, support, and understanding. Last but not least, I want to express my thanks to my 2-year-old son, Kevin, who always woke up at 2 am for milk, driving his sleepless father to focus on his work and thesis.   1  Chapter 1: Introduction  1.1 Overview of functional MRI In recent years, the technology of blood-oxygenation-level-dependent functional magnetic resonance imaging (BOLD-fMRI) has been highly developed. In addition to the advances of the devices and techniques of scanning, computer science related technologies such as digital images and graphics have brought about great advances in fMRI data analysis. The pre-processing and post-processing pipeline has become a critical component in fMRI data analysis.  A brief introduction of the fMRI data acquisition process and the data structure is given first, as shown in Figure 1-1. The patient for fMRI scan and analysis is called a subject. The fMRI scanner is a machine that could measure and map the brain’s activity of the subject. It’s a non-invasive test that uses a strong magnetic field and radio waves to create an image of the brain. In fMRI scanner, the brain of the subject is scanned slice by slice and finally we could get N 2-dimensional slices. The 2-D slices could also be called images as fMRI is an imaging technology. Then the slices are assembled one by one to create a 3-dimensional brain that is called a timepoint or a volume. The process that the fMRI scanner gets one timepoint of the subject could also be called a scan. As shown in Figure 1-2, a fMRI test is consisting of tens to hundreds of timepoints. Therefore, the raw fMRI data of a fMRI test is a 4-dimensional data including 3-dimentional volume data and the dimension of time. 2   Figure 1-1. Workflow of fMRI acquisition process. 3   Figure 1-2. Data structure of raw fMRI data.  During the fMRI acquisition process, a structural data is also captured by magnetic resonance imaging (MRI). The workflow of capturing the structural data is almost the same as that in Figure 1-1. The difference is that MRI emphasizes more on the structure of the brain, so it only has one timepoint with a much higher resolution than functional data. Figure 1-3 shows an example of one slice of the functional and structural data in the same subject.            (a) structural data (MRI)    (b) functional data (fMRI) Figure 1-3. Example of structural data and functional data. 4   1.2 Properties of fMRI data The structural volume is obtained by high-resolution T1, T2, and fast spoiled gradient echo (FSPGR) three-dimensional imaging. The processing of the functional volume is the key to fMRI data analysis. The instantaneous changes of cerebral cortical activity ask for a corresponding imaging sequence with two properties. First, the imaging sequence should be fast enough to record the cortical activity caused by a stimulation task. Second, it should be sensitive to the shortening effect of T2* produced by deoxyhemoglobin which is a product of cerebral oxygen metabolism. Imaging sequences like Echo Planar Imaging (EPI) and Fast Low Angle Shot (FLASH) can satisfy these conditions. As the signal to noise ratio (SNR) of EPI higher, it is widely used in sampling fMRI images nowadays.   Using a series of inverse gradients in frequency coding, EPI could generate all the signals required to construct an fMRI volume by a single stimulation. Based on the technology of Gradient Echo-Echo Planar Imaging (GRE-EPI) which is stimulated by a small angle, a series of images (from few to dozens) could be obtained in quite a short repetition time (TR). Those sampled images constitute a brain volume, and the time of each experiment block must be an integral multiple of TR when designing the experimental block of fMRI (Epoch/block Paradigm). Actually, the response of hemodynamics is a slow process. The signal has a downward period followed by a rise after the task stimulation, peaks in 4 to 8 seconds, then decreased slowly, and then recovers in 11 to 14 seconds. Currently, the experiments of fMRI mainly have three classes: Block Design, Event-Related Design, and Resting State. The properties of fMRI data are not the same under different designs. In general, researchers need to choose the design model for their own purposes. 5   Block experiment design divides tasks into blocks. Stimulation sequences consisting of the same tasks (such as navigation, control tasks, etc.) cause continuous activation of related brain areas. It will induce strong changes in BOLD signals and the changes will be collected by fMRI instruments. In general, the duration of one task block is 30 to 60 seconds, and the number of blocks in one fMRI scanning sequence is 6 to 12. The duration of a block should be appropriately chosen. If the length of the block is too long, the oxygen content in blood will not increase when it reaches its maximum, making the fMRI data unsatisfying. On the contrary, if the duration of the block is quite short, the oxygen content in the blood cannot return to the baseline, resulting in a decrease in the magnitude of BOLD signals.  In event related designs, each event is separated in time by an inter-stimulus interval ranging from few to 20 seconds. In this model, the stimulus tasks are discrete, and the duration between two stimulus is fixed. Event related designs attempt to measure transient changes in brain activity. Unlike blocked designs, stimuli are presented in a random order rather than an alternating pattern which offers a higher flexibility in experimental terms. The time interval should also be appropriately chosen. If the time interval is too long, it may cause distraction of the subjects and cause some changes in the brain. And if the time interval is quite short, it is difficult for the BOLD signal to return to the baseline. Generally, the time interval of event-related design should be controlled at about 12 seconds.  6  Resting-state fMRI (rs-fMRI) is used in brain mapping to evaluate regional interactions that happens in a resting or task-negative state, when an explicit task is not being performed. In this thesis, we mainly use the fMRI data under resting state.  1.3 Interference of fMRI EPI sequence could generate hundreds to thousands of images in a few minutes’ experiment with a very fast acquisition speed. The brain volumes at different time points constitute the time-series images of EPI. At the expense of image resolution, the size of a typical EPI image acquisition matrix is 64 by 64. Increasing the resolution of the acquisition matrix will prolong the sampling time and lead to more severe geometric distortions. Besides, EPI sequence images are quite sensitive to the external magnetic field environment. Thus, weak BOLD signals will have a large number of interference components.  1.3.1 Head movement during scanning Head motion is the main interference during fMRI scanning. Error! Reference source not found. and Figure 1-4 shows the volume-wised and slice-wised displacement caused by head motion.   (a) volume-wised displacement 7   (b) slice-wised displacement Figure 1-4. Displacement caused by head motion during fMRI scanning. Although various physical methods have been adopted, the effect of head movement is still difficult to eliminate entirely. Its side effects are far more than the mismatch in the fusion of the functional and structural images. The BOLD signal could be as weak as 0.5-2.0% when the magnetic field strength is 1.5 Tesla. While the signal difference between two adjacent voxels is usually larger than 10% and will even increase to 70% when getting close to the edge of the brain. Any tiny head motion will have a big effect on the analysis of brain activation and even cause a complete error to the detection of the brain activation. As a fMRI volume is a 3-D object, it could achieve a movement of six degrees of freedom which are [x, y, z, paw, pitch, roll], as shown in Figure 1-5(a). Therefore, the estimation of head movement of one timepoint could be separated as 3 translations and 3 rotations as shown in Figure 1-5(b). 8       (a) 6 degrees of freedom      (b) motion estimation Figure 1-5. Estimation of head motion.  1.3.2 Noise interference of devices and physiological activities Most of the interference caused by scanning equipment and physiological activities are high-frequency noise. Physiological activities like breathing and heartbeat will have a more significant impact on the detection of BOLD signals, especially when they are related to designed tasks. At the same time, as the change rate of the signal in the activation area is limited, spontaneous physiological activities will cause thermal noise and high time-frequency fluctuation. Besides, the instability of scanning devices will cause low-frequency drift.  1.3.3 Susceptibility and N/2 artifacts  The magnetic sensitivity of human brain tissues, especially near air-containing cavities like paranasal sinuses, are different. It will cause the inhomogeneity of local magnetic field which will lead to the geometric distortion of the reconstructed EPI image in the direction of phase encoding and make the functional area inaccurate.  0 50 100 150 200 250image-0.3-0.2-0.100.10.2mmx translationy translationz translation0 50 100 150 200 250image-0.200.20.40.60.8degreepitchrollyaw9    Figure 1-6. N/2 ghost artifact in fMRI scan.  Meanwhile, because of the inaccurate acquisition timing and the inhomogeneous static magnetic field, the echoes in k-space present a specific phase difference. EPI pulse sequences consist of a series of echoes. In zig-zag acquisition, every second echo is acquitted in an alternate direction. In the process of image reconstruction, even numbered echoes should be time-reversed in order to match the odd numbered echoes before Fourier transformation. If the forward and backward echoes are not perfect mirror images of each other, it will introduce artifacts. A simple delay of the start of the first echo will be propagated into all later echoes. It results in slight timing differences between the peaks of even and odd numbered echoes. Then when the Fourier transform 10  is performed, the phase error will result in signal intensity displaced in the phase-encoding direction halfway across the image. Assuming there are N pixels in the field of view, the aliased ghost appears shifted N/2 pixels relative to the main image positioned at the correct location. That’s where the name ‘N/2 artifacts’ comes from.   Figure 1-7. Conventional pre-processing pipeline for fMRI data.  1.4 Conventional pre-processing pipeline for fMRI data Targeting those interference in 1.3, a series of pre-processing procedures need to be taken to clean the raw fMRI data. Figure 1-7 shows a typical pre-processing pipeline for raw fMRI data. The 11  reason for taking each step and the objectives of those procedures will be detailed step by step in the next several sub-sections.  1.4.1 Realignment Realignment is always the first step of each pre-processing pipeline for fMRI data and it has two sub-procedures which are slice timing correction and motion correction. The acquisition time (tens of milliseconds) between EPI multi-layer acquisition in a single brain volume is slightly different. In block design experiments, as each task block lasts for a long time (from seconds to tens of seconds), these time differences can be ignored. However, in event-related design, the time requirement for task activation is relatively higher. It is necessary to correct the time differences collected from each layer to ensure that the tens of layers of one volume are completed at the same time. This operation is called slice timing correction. Head motion correction is a single modality-based registration of an ideal subject which is usually based on rigid body motion model (Grootoonk, 2000). Translation and rotation parameters are calculated iteratively to minimize the mismatch between the reference timepoint (usually the first timepoint of the time series) and the subsequent sequence timepoints, so as to implement the registration of all the time series timepoints. Three-dimensional space motion correction (Crinion, 2007) uses six parameters of translation in three directions and rotation in three coordinate axes to describe the head rigid body model. Hundreds to thousands of timepoints are collected in each experiment. A large amount of data makes the correction process quite time-consuming. In order to achieve real-time processing, some commercial software has abandoned this procedure. The development of fast motion correction algorithm is of great significance to real-time imaging.  12  1.4.2 Co-registration Low-resolution EPI functional images often need to be superimposed to high-resolution anatomical images to identify functional areas. It could be achieved by registering the activation mapping areas to the anatomical images. Unwrapping correction should be adopted for correcting the geometric and intensity distortion caused by ghost effect and magnetic sensitivity effect in this step.  1.4.3 Normalization The key of fMRI visualization is to map the detected functional activation areas to the high-resolution anatomical structure image accurately. The functional activation map does not contain any anatomical information and thus cannot be registered to the anatomical image. However, as the functional activation map and the functional image share the same coordinate system, it is possible to register the functional image and the anatomical image first, and the obtained transformation parameters can be applied to the functional activation map. The convolution parameters of the mean image generated by space correction/the registered anatomical image and the pre-designed template image in standard anatomical space are deployed to each slice. In this way, the image data of different subjects and modalities can be evaluated in the same coordinate system (Crinion, 2007).   13   Figure 1-8. Spatial normalization of fMRI volume.  For a single subject, the functional image could be normalized to an individual template rather than the standard space. For the brain structure highly damaged subjects like brain occupancy and infarction, a standard template should never be used. All the specific information of the damaged parts will be eliminated during the linear and non-linear transformation procedures. For these subjects, except for creating individual templates, cost-function masking could be used to cover the lesions first and then normalize the image to the standard space. Talairach and Tournoux system is the most classical standard anatomical system. The data comes from solid anatomy. The MNI system established by the Montreal Neurological Institute of McGill University was mapped 14  to Talairach and Tournoux using MR brain scans of 305 normal subjects. In this thesis, the MNI space is adopted.   Meanwhile, as the cerebral blood flow changes and the scanning hardware is not stable, the average signal intensity of the time series will have non-function-activity-related changes. It will make the response of each stimulation not at the same level and reduce the performance of analyzing the function activation information statistically. Each image needs to be adjusted so that their mean values are on the same level, which is the normalization of time series.  1.4.4 Smoothing The interference caused by hardware instability and physiological motion could be eliminated by smoothing as shown in Figure 1-9. Spatial smoothing can reduce the random noise of the MR image. It could also help to improve the signal-to-noise ratio and the ability to detect the functional activating data (Geissler, 2005). A filter could be generated by convoluting the fMRI data with a three-dimensional Gaussian function. The range of the smoothing filter could be expressed by the full width and half height (FWHM) of the Gaussian kernel.   In theory, the Gaussian kernel should have the same size as the reaction zone. At the same time, the Gaussian kernel must be larger than a voxel. Otherwise, it will cause data re-sampling and reduce internal resolution. When the signal-to-noise ratio is low, the detected activation region will cover a larger range if a wider filter is used. FWHM should also be larger when comparing multiple subjects so that the data of each subject can be projected to the same functional anatomy image 15  and reduce the difference between subjects. Although the filter can effectively remove the noise at a specific frequency, some real BOLD signal at the same frequency will be lost as well.   Figure 1-9. Spatial smoothing with Gaussian kernel.  1.4.5 Segmentation Brain segmentation is one of the most important tasks in medical image analysis and is often the first and the most critical step in many clinical applications. In brain MRI and fMRI analysis, image segmentation is commonly used for measuring and visualizing the brain’s structures and activation map. In this step, as shown in Figure 1-10, the brain is segmented into white matter (WM), grey matter (GM) and cerebrospinal fluid (CSF). The time series in different tissues are usually analyzed separately.  16   Figure 1-10. Brain segmentation.  1.4.6 Denoising The pre-processing procedures above are mainly spatial based. Denoising aims for cleaning the BOLD fMRI signal using data driven methods (time series) which are less specific in the assumed model and can simultaneously reduce multiple noise fluctuations. Most denoising approaches are based on data decomposition techniques such as principal and independent component analysis. In this thesis, we try to cover both conventional and novel denoising methods. 17  1.4.6.1 Realignment regression It has been a common strategy for alleviating the effects of head motion by considering the estimated motion traces in fMRI analyses. This assumes simple relations (linear or quadratic) of the residual motion-induced signal variations concerning detected displacements or motion velocity and appears to reduce the signal variance systematically (Lund, 2005). For both resting-state and task-based fMRI, motion parameters estimated in the realignment step are often used as additional confounds in the general linear model (GLM) analysis (Hutton, 2011; Friston, 1996). It has been widely used in some fMRI analysis toolboxes like SPM, FSL and AFNI to include the six realignment estimates, which is R = [x, y, z, pitch, yaw, roll] (6DoF), for regression since years ago. Later on, some subsequent studies have confirmed that the motion correction based on six parameters is inadequate. Therefore, for a better estimation, more motion-related covariates like the first and second temporal derivatives, [R, R'] (12 motion regressors), [R, 𝑅4, 𝑅567, 𝑅5674 ] (24 motion regressors where t-1 is the previous timepoint), and [R, 𝑅4, 𝑅567, 𝑅5674 ,	𝑅564, 𝑅5644 ] (36 motion regressors which uses the previous two timepoints), were introduced.  Typically, these approaches assume that head motion only occurs between volume acquisitions, implicitly ignoring the in-volume or slices misalignment. It just uses rigid body alignment and interpolation to eliminate the influence of detected motion. The differences between various realignment algorithms are that they use different interpolation algorithms and cost functions, as well as the optimization approaches to minimize (or maximize) the respective cost function.  1.4.6.2 Censoring For dealing with the incomplete removal of motion artifacts, another well-known method 'scrubbing' or 'censoring' has been used. In general, scrubbing excludes the individual high-motion 18  timepoints from the whole 4-D signal. Many groups have found that censoring improves the removal of motion artifact. Incorporation of global signal regression and censoring throughout processing could eliminate the functional connectivity correlations at nearly all distances (Power, 2012; Power, 2013; Power, 2014; Satterthwaite, 2013; Yan, 2013). However, some obvious drawbacks could be seen from this algorithm.  • Censoring(scrubbing) and spike regression will destroy the raw structure of the data, which will impact frequency filtering. • Direct removal of high-motion volumes will cause a high loss of temporal degrees of freedom(tDOF). • Scrubbing is vitally dependent on the estimation of the realignment parameters derived from the motion corrected fMRI volumes, which may not be accurate and therefore cannot be regarded as a baseline for scrubbing.  Another concern of censoring is that it is difficult to decide the threshold of motion artifacts. In order to censor more motion-contaminated data, a lower threshold needs to be set. Thus, it introduces a balance game between censoring data to reduce motion-related artifact and the demand of preserving the data as much as possible.  1.4.6.3 aCompCor Another popular denoising method is aCompCor (Behzadi, 2015). It defines areas whose signal is presumably of no interest and applies temporal principal component analysis to derive nuisance regressors. By extracting orthogonal components of time series fMRI signal from WM and CSF, 19  aCompCor has a better ability to characterize the noise components in those brain tissues. In our comparison experiments, five principal components are extracted in aCompCor.  1.4.6.4 ICA-AROMA Independent component analysis (ICA) based denoising approaches have been widely applied. As a blind source separation (BSS) technique, ICA could identify maximally spatially independent components in the signal in theory. Recently, an ICA based motion denoising method, ICA-AROMA was introduced (Pruim, 2015; Pruim, 2015). This approach could identify and remove motion-related artifacts from raw fMRI signal automatically based on matrix decomposition. Independent components are classified as motion related or as useful signal according to several criteria.  ICA-AROMA is a subject based denoising method as well and needs to be applied to each scan separately. As the thresholds of the classifier for detecting motion-related components are set based on a specific subject space, the number of the confound regressors are always not the same across the whole dataset. Besides, the spatial smoothing should be applied before using ICA-AROMA. Otherwise, the signal of fMRI voxels could be discrete and nearly all the components derived from ICA could be regarded as motion-related and be removed.  1.5 Effect of head motion Head motion has a considerable effect on the analysis of fMRI signals while analyzing both task- related and resting-state fMRI data. Motion artifacts could still contaminate the raw BOLD data even if they are quite small when compared to the voxel resolution (e.g., 3 × 3 × 3	mm;). It makes 20  it even harder to decide whether the intensity variation of the BOLD signal comes from real brain activities or is just due to head motion. Eliminating the influences of head motion has been an essential issue in fMRI data analysis since two decades ago (Jenkinson, 2002).   Most motion correction methods use a rigid-body setting (assuming that the brain volume is constant between frames) and only correct for translation and rotation. Typically, an implicit assumption is that head motion only occurs all at once between two volume acquisitions, ignoring the inter-volume or slices displacement when motion occurs. However, in fact, considering how motion displacement actually is generated, the effect on the BOLD signals cannot be merely linear.  If such head motion effects are not sufficiently eliminated by preprocessing, the remained spurious changes can significantly confound the estimation of the functional connectivity between brain regions learned from the fMRI data.               Figure 1-11. Illustrative examples of correlation changes caused by head movement.  Red connections indicate the increased correlations and blue connections indicate the decreased correlations between ROIs, when compared with the corresponding correlations using the original, clean fMRI data.   21  Unlike in task-based fMRI signals, motion artifacts in resting - state fMRI could not be suppressed by averaging. Depending on the non-neural signal fluctuations that subject’s motion generates in different ROIs, motion-related fluctuations can cause both increases and decreases in the correlations between those ROI fMRI time series. Therefore, head movement can create spurious correlations at both the subject level and the group level, and it portends something more challenging for the studies of functional connectivity (Power, 2014; Van Dijk, 2012). When head motion occurs, two ROIs may experience distance-dependent modulation of fMRI signal correlations. More specifically,  such a modulation can 1) increase the correlations between nearby ROIs; 2) increase the correlations between some distant ROIs when the motion is widespread across many ROIs; 3) decrease the correlations between distant ROIs (occasionally also nearby ROIs) when the motion shows a different property at different locations of the brain. (Power, 2015) (Power, 2012). Accordingly, the analysis and results of functional connectivity using fMRI may become unreliable. In some clinical studies, this situation is exacerbated when dealing with cases that the scanning time is limited and repeating scans are not always possible, for example, for patients with Parkinson's disease (PD) or children and infants.  1.6 Motivation, objectives and thesis organization In resting-state fMRI, though there is no stimulation of externally prompted task, the brain regions will have spontaneous fluctuations in BOLD signal because of the intrinsic brain activity. The resting-state fMRI is extremely helpful to examine if the functional organizations of the brain are 22  altered in neurological or mental disorders. The research of resting-state functional connectivity has also revealed many networks which are consistently found in healthy subjects, different stages of consciousness and across species, and represent specific patterns of synchronous activity (Rosazza, 2011).   However, as resting-state fMRI evaluates regional interactions that occur in a resting or task-negative state, the BOLD signal has lower signal to noise ratio compared with task-related fMRI. Therefore, the BOLD signal is much more sensitive to the distortion of the brain volume caused by head motion. An effective motion correction method is significantly important to the study of resting-state fMRI.   Besides, we believe that if the proposed motion correction method could help to remove the motion artifacts in resting-state fMRI, it could be easily applied to task-related fMRI as well. That is the reason why we focus on the resting-state fMRI data.   According to the introduction in section 1.3 and 1.4, conventional motion correction methods have some deficiencies: 1. The result of motion correction is highly dependent on the accuracy of motion estimation,which increases the accumulative error of the pipeline. Actually, in fMRI analysis, it’s impossible to get an accurate motion estimation result because any tiny change of the brain including both head motion and physiological activities will change the result of motion estimation.   23  2. Conventional motion correction methods treat the brain as a rigid body, ignoring the inter-volume or slice-wised displacement. However, during the acquisition of fMRI, the time interval between two timepoints only takes a small fraction of the whole process. It means that inter-volume displacement occupies the main part of the motion artifacts.  3. During the rigid-body based motion correction, the interpolation is always used. But no matter what the strategy of the interpolation is (linear, cubic or spline), the intensity value of the voxel is created based on its surrounding voxels. While the signal difference between two adjacent voxels is usually larger than 10% and will even increase to 70% when getting close to the edge of the brain. It is a compromise on calculation rather than the estimation of the real voxel.  To tackle the shortages of conventional motion correction methods of fMRI data, we aim to develop a novel IVA-based method for efficiently removing artifacts caused by head motion.  The organization of this thesis is as follows:  Chapter 2: Motion correction with IVA. In this chapter, the independent vector analysis concept is introduced, and the proposed motion correction model of McIVA is detailed.  Chapter 3: Evaluation of McIVA. In this chapter, we mainly evaluate the performance of the proposed McIVA model on both motion correction and denoising. The dataset and preprocessing methods are described first, and several assessment measures are introduced.  24   Chapter 4: Results of experiments. In this chapter, the assessment results are reported.  Chapter 5: Discussion.  In this chapter, some further analysis of the proposed algorithm and the structure of McIVA, as well as the results in chapter 3, are discussed.   Chapter 6: Conclusion and future work.  In this chapter, we summarize the contributions of this thesis and discuss the future direction of this research.25  Chapter 2: Motion correction with IVA  To address current issues associated with non-linear displacement of head movement and rigid body correction, we propose an alternative IVA based strategy which is called 'McIVA' or ' Motion Correction with Independent Vector Analysis'.  2.1 Motivation The target of motion correction is to remove the motion related effects with the prerequisite of not eliminating the active information of the brain. To a 4-D fMRI signal, the correlation of all the volumes should be the highest when there is no motion occurred. As a practical matter, during the test of one subject, head movements will affect all the scans more or less. However, the information shared across all the scans of the subject and all the sub-volumes of the motion corrupted volume is the essential volume information of the brain with no motion artifacts. Unlike ICA-based method whose target is to remove the most motion-related signals, the goal of McIVA is to extract the essential feature or the non-motion-contaminated volume from the sub-volumes under the principle of maximizing the correlation to the reference volumes.   Meanwhile, as all the sub-volumes are generated from the motion corrupted volume, it will preserve the useful information of the brain. The volume with a red frame in the fMRI time series is motion contaminated volume. The module in the red dashed frame is the IVA core to extract the motion-corrected volume. Observation volumes are generated through rigid body transformation of the motion corrupted volume. The number of observation volumes M should always be larger than the number of reference volumes N. 26  2.2 Independent Vector Analysis Independent vector analysis (IVA) is a multivariate extension of independent component analysis (ICA), addressing the limitation of the conventional ICA method during blind source/signal separation (BSS) (Kim, 2006).   BSS could be conducted in both time domain and frequency domain, and due to the utilization of FFT, BSS on frequency domain could offer a faster convergence speed and less computational demand. However, in the frequency domain, BSS may result in the permutation of the identified IC arrays (within each frequency bin), which requires additional steps (such as correlation analysis) to search for the mutually dependent components across frequency bins. IVA was proposed to tackle this permutation problem in frequency-based BSS to correctly index these independent components. Instead of a scalar value in ICA, each component of the input and output stages in IVA forms a vector. IVA tends to increase independence among those vectors under the restriction of keeping the dependency of the elements within each IC vector.  The mathematical expression of the linear mixed model of ICA could be expressed as, 𝑥 = 𝐴𝑠 = 𝐷𝑃𝑠 (2.1) where s = (𝑠7, 𝑠4, … , 𝑠I)J  is N independent sources, 𝑥 = 	 (𝑥7, 𝑥4, … , 𝑥I)J  is the linear combination of s. A is a constant mixed matrix of N × N, (∗)J is matrix transposition, and both s and A  are unknown. D  is a diagonal matrix indicating the amplitude ambiguity, and P  is a permutation matrix indicating the permutation ambiguity. The target of the ICA algorithm is to estimate the demixing matrix 𝑊, and use the demixing matrix to estimate the source signal y, where 27  y = W𝑥 (2.2) While IVA combine multiple 𝑥 and make a separation to all those 𝑥, which could be expressed as 𝑥(Q) = 	𝐴(Q)𝑠(Q),			m	 = 	1, 2, 3, … ,M (2.3) in which, M is the number of the array, 𝑠(Q) = [𝑠7(Q), 𝑠4(Q), … , 𝑠I(Q)]J and 𝐴(Q) are the sources and mixing matrix of mth array. The nth source component vector is defined as 𝑠S =	[𝑠S(7), 𝑠S(4), … 𝑠S(T)]J.  2.3 IVA-L algorithm IVA-L algorithm sets the joint distribution of the original signal as a multi-dimensional Laplace distribution and uses divergence as the cost function. The source component vectors are constrained to be independent by minimizing the divergence. IVA-L estimates the demixing matrix 𝑊(Q),𝑚 = 1,… ,𝑀, and uses the demixing matrix to calculate the sources 𝑦(Q) = 𝑊(Q)𝑥(Q), and then the nth source component vector could be derived which is 𝑦S = [𝑦S(7), … , 𝑦S(T)]J. The cost function of IVA-L is shown as below. 𝐼 = 𝐾𝐿 [𝑝(𝑦7, … 𝑦I) ∥^𝑝(𝑦S)IS_7 ` = 	a𝑝(𝑦7, … 𝑦I)𝑙𝑜𝑔 𝑝(𝑦7, … 𝑦I)∏ 𝑝(𝑦S)IS_7 𝑑𝑦7 …𝑑𝑦I (2.4) where 𝑝(𝑦7, … , 𝑦I) is the joint density function, and p(𝑦S) is the marginal density function of single SCV. When the distribution of SCV is independent, the joint density function of the SCVs is the product of the marginal density function, and the minimum value of I is 0. The optimization method of IVA-L is natural gradient algorithm and the iterative function could be expressed as 𝑊(Q) 	→ 	𝑊(Q) + 𝜂Δ𝑊(Q) (2.5) 28  where  𝑊(Q) is the demixing matrix of mth array, 𝜂 is the iteration parameter. The model IVA-L used is multidimensional Laplacian probability density model, whose distribution could be expressed as p(y) = Cn 𝜋2√2𝜆 exp s−u2𝜆(𝑦 − 𝑢)JΣ67(𝑦 − 𝑢))x (2.6) Where C is a constant, u is the mean value of vector y, Σ is the covariance matrix of y. The prior condition in IVA-L is that the source vector is an independent vector with zero mean, therefore, u = 0, Σ = I. Then we have p(y) = Cexp}− ~|𝑦|4TQ_7 € (2.7) According to equation (2.4) and using the natural gradient algorithm, ∆𝑊(Q) could be derived as ∆𝑊(Q) → [𝐼 − Φ(𝑦(Q)(𝑦(Q))J)]𝑊(Q) (2.8) in which Φ(∗) is the score function,  Φ…𝑦(Q)† = ‡Φs𝑦7(Q)x, … ,Φs𝑦I(Q)xˆJ (2.9)  moreover, Φs𝑦S(Q)x = − Š‹ŒŽ()ŠŽ((‘)) . After putting equation (2.7) into (2.9), we could get Φ…𝑦(Q)† = ’−𝜕𝑙𝑜𝑔𝑝(𝑦7)𝜕𝑝s𝑦7(Q)x , … ,−𝜕𝑙𝑜𝑔𝑝(𝑦I)𝜕𝑝s𝑦I(Q)x ” (2.10) According to equation (2.8) and (2.10) the demixing matrix 𝑊(Q)  of the mth array could be derived.   29  2.4 IVA-G algorithm There are several kinds of IVA-G algorithm according to different optimization algorithms. The multidimensional density model of IVA-G-M is multidimensional Gaussian distribution, and the optimization method is matrix gradient. The cost function uses mutual information and could be expressed as 𝐼 ≜ 𝐼[𝑦7, … , 𝑦I]J = 	~𝐻[𝑦I] − 𝐻[𝑦7, … , 𝑦I]IS_7 (2.11) where 𝐻[𝑦I] is the entropy of 𝑦I . According to the property of linear conversion of entropy, H[Wx] = log	 |𝑑𝑒𝑡(𝑊)| − 𝐻[𝑥], where det	(∗) is the determinant of a matrix. And according to the definition of mutual information, H[𝑦S] = 	∑ 𝐻‡𝑦S(Q)ˆ − 𝐼[𝑦S]T . Then we have the cost function of IVA-G-M which could be expressed as I = 	~}~ 𝐻‡𝑦S(Q)ˆ − 𝐼[𝑦S]TQ_7 €IS_7 − ~ logdet…𝑊(Q)†TQ_7 − 𝐶𝑜𝑛𝑠𝑡 (2.12) Where N is the number of the SCVs, and M is the number of the observation datasets. According to equation (2.12), we could see that by minimizing the cost function, the entropy of the source vectors could be minimized and meanwhile the mutual information among the SCV could be maximized.  Multidimensional Gaussian model was used as the joint probability density distribution in IVA-G-M which could be expressed as p(𝑦S|ΣS) = 	 1 2𝜋T4 ¡ det(ΣS)74 exp  −12𝑦SJ~ 𝑦S67S ¡ (2.13) 30  where ΣS is the covariance matrix of second order statistical information of SCV. As the second order statistical information was used, the problem of ambiguity could be eliminated better. Then we take the derivative of equation (2.12) Δ𝑊(Q) = 	−~𝐸 £𝜕log𝑝(𝑦S)𝜕𝑦S(Q) 𝜕𝑦S(Q)𝜕𝑊(Q)¤IS_7 − …𝑊(Q)67†J (2.14) Then according to equation (2.9), we have Φ(𝑦S) = 	∑ 𝑦S67S .  Take this equation to (2.14) to get Δ𝑊(Q) = 	𝐸 ¥ΣS67𝑦S…𝑥(Q)†J¦ − …𝑊(Q)†6J (2.15) According to the natural gradient, the iteration function could be gotten 𝑊(Q) ← 𝑊(Q) − 𝜇 s𝐸 ¥ΣS67𝑦S…𝑥(Q)†J¦ − 𝐼x𝑊(Q) (2.16) The IVA algorithm IVA-L utilizes a multivariate circular Laplace prior as a distribution on the SCVs. Thus, the components within each SCV are assumed to be second-order uncorrelated and IVA-L exploits only higher-order statistics. On the contrary, IVA-G exploits only second-order statistical information by multivariate Gaussian assumption. The Matlab implementations of these algorithms are available at http://mlsp.umbc.edu/resources.html.  In this thesis, the IVA-G algorithm was adopted in McIVA to extract the first correlated independent vector from the subvolumes of the motion corrupted volume. The choice of IVA algorithm didn’t make a big difference in the performance of McIVA.  2.5 Model structure The main structure of McIVA is shown in Figure 2-1. The volume with a red frame in the fMRI time series is motion contaminated volume. The module in the red dashed frame is the IVA core 31  to extract the motion-corrected volume. Observation volumes are generated through rigid body transformation of the motion corrupted volume. The number of observation volumes M should always be larger than the number of reference volumes N.   Figure 2-1. The main structure of the proposed McIVA. The volume with a red frame in the fMRI time series is motion contaminated volume. The module in the red dashed frame is the IVA core to extract the motion-corrected volume. Basis volumes are generated through rigid body transformations of the motion corrupted volume. The number of basis volumes M should always be larger than the number of reference volumes N.  Some magic numbers are predefined to have a better demonstration. The 4-D fMRI signal which is consist of 3-D fMRI information and the time dimension is one test of No. K subject. Each time point in the test is called a scan or a volume. N is the number of the reference volumes or called the sources volumes. M is the number of sub-volumes of the motion corrupted volume or called the observation volumes. V is the number of voxels in one scan.  32   N volumes are first chosen as the reference volume based on the principle of minimal proximity to the mean volume. Then each of the reference volumes is transformed to a 1	 × 	𝑉 vector to make the source matrix 𝑠I׫	 = [𝑠7	, 𝑠4	, …	𝑠I	],  where 𝑠I	 = 	 𝑠I(1), 𝑠I(2), … 𝑠I(𝑉). In fact, as the effect of head movement spreads in the time dimension, the reference volumes, though chosen carefully, still contain motion artifacts. So McIVA is first applied to each of the reference volumes, taking one as 'motion corrupted volume' and the other N-1 volumes as reference volumes.  For the motion contaminated volume, M sub-volumes are extracted first through rigid body transformation. Then each of the sub-volumes is transformed to a 1	 × 	𝑉 vector to make the initial observation matrix ?⃗?T׫	 = [?⃗?7	, ?⃗?4	, …	?⃗?T	], where  ?⃗?T	 = 	 𝑥T(1), 𝑥T(2), … 𝑥T(𝑉). The number of the sub-volumes of the motion contaminated volume M is expected to be large enough to contain more features under different degrees of translation and rotation. Subsequently, the principle component analysis is applied to extract the top N principle components to build the final IVA observation matrix 𝑥′­­­⃗ I׫	 = [𝑥′­­­⃗ 7	, 𝑥′­­­⃗ 4	, …	𝑥′­­­⃗ I	]. The core IVA algorithm is then applied to the two matrices 𝑠I׫	 and 𝑥′­­­⃗ I׫	.   The advantage in the ordering of the IVA algorithm makes it possible to extract the principle vector which is the most relevant to the reference vectors. According to the analysis we made at the beginning of this chapter, the property of correlation is a necessary but insufficient condition of the effectivity of motion correction. More specifically, the correlation of one volume and the reference volumes reaches a maximum value when it contains no motion artifacts. Nevertheless, 33  it is imprecise and inaccurate to conclude that when the correlation reaches the peak, the corresponding volume is non-motion corrupted. Intuitively, a volume with brain activation has a lower correlation value than the one with no activation. While in McIVA, an implicit domain restriction is applied to the volume to be corrected, which is that the extracted volume comes from the space consisted of the sub-volumes of the motion corrupted volume. Thus, a conclusion could be drawn that within the space of the sub-volumes of the motion contaminated volume, the extracted volume which maximizes the correlation to the reference volumes has the least motion artifacts.  For higher accuracy, an adaptable method is adopted according to the degree of the estimated motion when extracting the sub-volumes from the motion corrupted volume. In general, the step of the rigid body transformation is set relatively larger when the estimated motion is bigger and vice versa. The purpose of taking this step is to acquire as much space features of the motion corrupted volume as possible. In this way, the volume with no motion effects could be extracted from the space consisted of the sub-volumes.     34  Chapter 3: Evaluation of McIVA  3.1 Dataset and Preprocessing 3.1.1 Participants Thirty-two subjects with PD (10 females; age: 66.22±6.62; UPDRS III: 23.56±12.00; Hoehn and Yahr scale: 1.96±0.93) and fifteen age-matched healthy controls (HC, five females; age: 69.4±4.76) were recruited from Pacific Parkinson’s Research Centre (PPRC) at the University of British Columbia (UBC). The PD subjects had mild to moderate PD (Hoehn and Yahr stage I - III) and were clinically in the on-medication state. The study was approved by the UBC Ethics Review Board, and all the participants provided written, informed consent before the experiment.  Table 3-1. Characteristics of the participants included in the dataset.  N AGE (years; mean±	SD) Female (%) PD 32 66.22 ± 6.62 31.25 HC 15 69.4 ± 4.76 33.33   3.1.2 Data acquisition A 3T scanner (Philips Achieva 3.0T R3.2; Philips Medical Systems, Netherlands) equipped with a head-coil was used to obtain the resting-state imaging data. Before scanning, all the subjects were instructed to lie on their back in the scanner and have several minutes to acclimatize themselves to the scanner environment with eyes closed. High-resolution T1 weighted anatomical images were acquired with a repetition time of 7.9 ms, echo time of 3.5 ms, and a flip angle of 8°. 35  Blood oxygenation level-dependent (BOLD) contrast echo-planar (EPI) T2*-weighted images were acquired using with the following specifications: repetition time of 1985ms, echo time of 37 ms, flip angle of 90°, field of view of 240 mm × 240 mm, matrix size of 128 × 128, pixel size of 1.9 mm × 1.9 mm and the scanning time of 8 mins for rest condition and 5 mins for GVS conditions, respectively.  3.1.3 Data preprocessing For the preprocessing of the Brain fMRI data, the CONN toolbox (https://www.nitrc.org/projects/conn) (Whitfield-Gabrieli, 2012) SPM12 toolbox (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) (Penny, 2011) and FSL (https://fsl.fmrib.ox.ac.uk/fsl/fslwiki) (M.W. Woolrich, 2009) were used.   The dataset was preprocessed in two workflows to form two datasets for different purposes. For dataset-1, only slice time correction and realignment were applied for comparing the performance of McIVA with common motion correction methods.  For dataset-2, the dataset was processed as the workflow shown in Figure 3-1 on each subject. First, slice time correction was applied to the fMRI data to correct for different scan timings and shift the time courses so that all slices appear to be acquired simultaneously. Second, for non-McIVA methods, functional volumes were realigned to a subject-based global mean for motion correction via a rigid body transformation. For McIVA, the McIVA method and motion estimation were taken. The alignment and estimated parameters were saved as confounds and will be regressed out for nuisances regression approaches. Third, the center of the functional volumes was  36  translated to (0, 0, 0) coordinates. Forth, ART-based identification of outlier scans was applied for scrubbing. Fifth, simultaneous Gray/White/CSF segmentation was applied to all the subjects, and the fMRI volumes were normalized to MNI space. Finally, a spatial convolution with a Gaussian kernel was applied.   Figure 3-1. Preprocessing pipeline for dataset in Table 3-1.  37  For the structural volumes, first, the centers were translated to (0, 0, 0) coordinates. Second, simultaneous Gray/White/CSF segmentation and normalization to MNI space are applied to all the MRI volumes.   3.2 Measure methods 3.2.1 Absolute frame displacement In order to test the performance of motion correction on the subject level, initial tests used the absolute displacement which indexed the absolute displacement of the head from the mean position at every time point in combination with volume based mutual information. The absolute displacement had six parameters: translational displacement of X, Y, Z and rotational displacement of pitch, yaw, and roll. Combined with the translational and rotational displacements, a QC measure was taken which converted the rotational displacements to translational displacements by projection to the surface of a 50mm radius sphere. The absolute displacement was the summation of the six absolute values.   3.2.2 Mutual information When revising the influence of motion artifacts, it should be ensured that motion correction does not remove the useful raw information of rs-fMRI. As a robust method in quantifying the relationship between two waveforms, mutual information (MI) has been used for analyzing fMRI since years ago (Tsai, 1999). From the perspective of information theory, the increase of mutual information indicates an increase of the correlations among the volumes in the test. Here the volume based mutual information was evaluated to estimate whether mutual information between volumes will be improved after motion correction. 38   3.2.3 Assessment of interpolation effects Some recent research found that processing procedures including temporal interpolation will have a significant effect on motion estimation. By applying temporal interpolation, estimated head movement was reduced by 10-50% or even more which makes the data to have a falsely improved quality (Power, 2017). Many motion correction methods are implicitly or explicitly based on a two-step process. The first step is to estimate the motion accurately. Then the motion is corrected by an interpolation process to reconstruct the image in its presumed original position. While much effort has been put into accurate motion estimation schemes, errors induced by interpolation are generally assumed to be inevitable, and hence ignored. However, since the proposed scheme takes a linear combination of basis volumes (each volume itself subject to interpolation errors), it is possible that the final, registered volume will have a different magnitude of errors compared the errors introduced by interpolation.  A volume from a real fMRI time series was selected and initially copied to ensure no motion. We then systematically rotated and translated in the X, Y, Z directions and around these axes, with parameters [2:0.5:2] mm and applied cubic interpolation with the Matlab command ‘maketform’. Since we knew the transformation, we were able to apply the exact inverse transformation subsequently. If interpolation error were zero, the mean absolute error in voxel intensity between the original motionless image and the transformed/inverse transformed image would be zero. We also corrected the transformed images using McIVA for comparison. The experiment of assessing interpolation effects was carried out on dataset-2. One complete brain volume with less distortion was chosen and extracted manually from one subject as the original 'no motion' volume. 39   3.2.4 Outlier detection Sometimes the data are corrupted that it is virtually impossible to reconstruct a motion-compensated image. We therefore purposely created an extremely corrupted volume to see if there was a rigorous way to detect such outliers. We took a middle volume in the series, (the 50th one). It had an intensity range from 0 - 1300. Then the intensity of the voxels was rearranged randomly under the premise of keeping the global mean value and the grey histogram unchanged, as shown in Figure 3-2. The experiments for outlier detection was carried out on one subject in dataset-2. The 50th volume was extracted when the head motion was relatively small.             (a) Original volume   (b) Artificially corrupted volume Figure 3-2. The timepoint extracted from real fMRI data and artificially corrupted volume.  3.2.5 Functional connectivity correlations The Pearson correlations are applied to quantify the association of subject-specific displacement and the estimation of functional connectivity, resulting in a QC-FC estimation on both subject level and dataset level. Head movement could not only decrease the correlation among volumes 40  but have an impact on the correlation among the ROIs within one volume as well. In addition, this impact is non-linear and shows a property of distance based partial correlation, which lies in that the degree of motion effect is not the same at different distances between ROIs. Some researchers have demonstrated that the motion will cause an entire increase in the functional correlation (Van Dijk, 2012). Thus, an critical evaluation of the performance of different motion correction approaches is to identify if they could decrease the functional correlation of those ROIs in the motion contaminated scans. In our experiments, both voxel-to-voxel and ROI-to-ROI based functional connectivity correlations were calculated for different purposes. Using dataset-1, we performed different motion correction schemes and then extracted the mean time courses from 192 ROIs before assessing the ROI-to-ROI correlations. We then chose a ‘seed’ region, in this case, the right insular cortex, and examined the correlation between this ROI and all other ROIs across subjects. This region was chosen as one where widespread connectivity would not necessarily be expected based on biological considerations (Poldrack, 2007).  We also explored the impact of motion correction approaches on different denoising schemes using the toolbox CONN. We first applied McIVA and then regressed out the mean signals from both CSF and WM before assessing voxel-to-voxel correlations across all volumes. This was contrasted to aCompCor, where SPM motion correction was used, followed by regressing out 5 Principal Components from WM and CSF as well as the 12 Realignment parameters. Finally, we also compared the results to IC-AROMA after performing SPM motion correction.  A bandpass filter of 0.008Hz to 0.09Hz was applied during the denoising procedure of the three approaches. For different purposes, the evaluation for functional correlation was separated into 41  three parts. The experiment for ROI-based correlation was carried on one subject in dataset-1 first to eliminate the interference of other preprocessing procedures. And then the experiments for both ROI-based and voxel-based correlation were carried out on dataset-2 to evaluate the performance of the whole preprocessing procedures. 42  Chapter 4: Results  4.1 Mutual information between first and subsequent volumes As anticipated, the MI between the first and subsequent volumes in the non-motion corrected data slowly decreased over time (Figure 4-1, blue line). After SPM motion correction (Figure 4-1, red line) the MI decreased approximately the same as the uncorrected data until time point 120 when it actually decreased faster. In contrast, after McIVA correction (Figure 4-1, green line) the mutual information decreased at a much slower rate. This experiment was carried out on dataset-2.  Figure 4-1. Mutual information between the first and subsequent volumes.  0 20 40 60 80 100 120 140 160 180 200 220timepoints0.720.740.760.780.80.820.840.860.880.90.92mutual informationRaw dataSPM correctedMcIVA corrected43  The trend of the mutual information shown in Figure 4-1 actually accords with the prevailing situation when sampling the fMRI signal. At the beginning of the process, it is easy for the patient (or subject) to stay dormant and relaxed. The mutual information between the first several volumes is higher than that of the following volumes. With the process going on, it becomes harder and harder for the subject to keep static especially for the patients with PD or children. The subject begins to get impatient. The mutual information decreases rapidly as head motion increases. Meanwhile, the brain also becomes more active which led to a further decrease in the MI value.   4.2 Mean mutual information between two volumes This part details the comparative accuracy of McIVA when tested against two of the most widely used toolboxes, SPM12 (released October 2014) and FSL6.0 (released October 2018). For SPM12, the default six realignment parameters were used for motion correction and for FSL6, a specific 6 and 12 parameters were adopted. This experiment was carried out on dataset-1. The reason that led to different results of SPM and FSL was that they took different cost functions for motion correction.  The MI between every pair of volumes in one subject was calculated to estimate the mutual information on a larger scale. As shown in Figure 4-2, all the four methods helped to improve the mutual information of the 4-D fMRI signal. Among the four methods, McIVA gives the best performance on improving the MI with a higher mean mutual information and lower variance. This is predictable as, unlike other methods, McIVA takes the conditional correlation-based optimization as its cost function, which gives it an advantage on reducing the otherness caused by motion artifacts of the volumes. 44   Figure 4-2. Mutual information of the original and the motion corrected fMRI data. The motion correction methods used in this figure is SPM, Dof6 and Dof12 of FSL and McIVA.  The difference in Figure 4-2 is not that obvious as shown in Figure 4-1. Figure 4-1 was the mutual information which changed with the time while Figure 4-2 was the MI of random two volumes which estimated the mean MI on subject level. Two volumes that both have a large displacement could have higher MI as well.   4.3 Mean displacement The mean value and the standard deviation of the absolute displacement on the subject level were calculated and demonstrated in Figure 4-3. Some group also use framewise displacement (FD) as a parameter which measures the displacement of every two adjacent volumes. However, an unusual situation is that, when the motion is a relatively long-term effect which lasts from the beginning to the end of the process, the FD could be extremely low, but the absolute displacement 45  of two non-adjacent volumes could be quite high. Therefore, the absolute displacement to the mean volume was calculated here to evaluate the motion correction performance. This experiment was carried out on dataset-1.  For absolute displacement to the mean volume, SPM12 shows a small improvement from 1.4mm to 1.18mm and two FSL6 methods show a raise after applying the motion correction. Among the four methods, McIVA shows the most significant improvement which decreases the absolute displacement from 1.4mm to less than 0.8 mm.   Figure 4-3. The average displacement of the original and the motion corrected fMRI data. The subject-level average displacement of the original fMRI data and the motion corrected fMRI data using different strategies. The motion correction methods used in this figure is SPM, Dof6 and Dof12 of FSL and McIVA.  46  4.4 Interpolation Error The effects of interpolation error were shown in Figure 4-4. After rotation and translation of a volume by varying degrees, there was some mean absolute error (MAE) between the original volume and the transformed volume (dotted black line). When the exact inverse transformation was applied to the transformed volume, as expected, the MAE was reduced (red line), but not to zero, due to interpolation errors.  After McIVA, the MAE was also reduced (blue line) and less than that of the interpolation error in many cases.  Figure 4-4. Mean absolute error of two approaches. The mean absolute error of the fMRI data motion corrected with the proposed McIVA and the traditional interpolation.  0 20 40 60 80 100 120 140Transformations12131415161718192021Mean absolute errorInterpolationMcIVAEstimated motion47  When examined at the individual voxel level (one volume was shown in Figure 4-5), the error with McIVA (blue dots) seemed to more consistent, compared to the interpolation errors, which varied substantially across the volume (red dots).   Figure 4-5. The voxel intensity changes of two approaches. The voxel intensity changes of the fMRI data motion corrected with McIVA and the traditional interpolation.  4.5 Evaluation of processing mechanism of outliers After severe corruption of the 50th volume SPM still "corrected" this volume (Figure 4-6). However, the IVA step with McIVA failed to converge for this particular time point suggesting that no matter how the corrupted volume is spatially transformed, it could not "match" the rest of the time series (estimated by the reference volumes). According to the theory of independent vector analysis, the independent vector which is most correlated to the reference volumes should be extracted. When the 50th volume was corrupted, IVA core was not able to extract an independent 48  vector with the corresponding cross-correlation coefficient. The forced extraction of the independent component will exacerbate the motion estimation of the corrupted volume, as shown in Figure 4-7.   (a) Original       (b)SPM Figure 4-6. Motion estimation of the original and motion corrected fMRI data with SPM. The estimated motion parameters were separated as translations (x y z) and rotations (pitch roll yaw). (a) Motion estimation of raw fMRI data with 50th volume damaged; (b) Motion correction  0 50 100 150 200 250image-0.3-0.2-0.100.10.2mmx translationy translationz translation0 50 100 150 200 250image-0.200.20.40.60.8degreepitchrollyaw0 50 100 150 200 250image-0.02-0.015-0.01-0.00500.0050.010.0150.02mmx translationy translationz translation0 50 100 150 200 250image-0.04-0.03-0.02-0.0100.010.020.030.04degreepitchrollyaw49   Figure 4-7. Motion estimation of the fMRI data motion corrected with McIVA.  Figure 4-8. Cross-correlation coefficient in McIVA. It shows the cross-correlation coefficient of reference volumes and independent vectors in McIVA.  0 50 100 150 200 250Index of volumes0.40.60.81Cross-correlation coefficient50  4.6 Functional connectivity In Figure 4-9, the ROI-to-ROI correlation was calculated after performing the four motion correction methods separately. Compared with the original correlation drawn in black, McIVA had a comparable better performance than the other three methods and was the only method that could decrease the ROI based functional correlation from middle distance to long distance. When two ROIs are close to each other, they will suffer a quite similar distortion caused by motion. Thus, the correlation between them will ascend rapidly even under a slight disturbance, which makes it extremely difficult to eliminate the motion effect in this situation. That is why the proposed method seems less effective when the ROI distance is small. Actually, it is difficult to use a single technique to remove all the motion effects and establish that the fMRI signal is free of motion artifacts. The advantage of McIVA lies in the fact that it could mitigate the difference on the basis of maintaining the inherent characteristics of the volumes. This technique is attractive because it does not rely on the features of the whole dataset and thus it could be applied to single scans. The benefits of its superiority will be revealed gradually in subsequent processing. 51   Figure 4-9. ROI based function correlation under specific motion correction methods. Head movement will have a local or global effect on the whole brain, leading to a positive distance-related shift of the correlation of the regions of interest. The unit of the horizontal axis is the number of voxels (3mm×3mm×3mm), and the vertical axis is the ROI-FC range from -1 to 1.  The advantage of McIVA lies in the fact that, although McIVA is applied based on rigid body estimation first, its degree of accuracy does not rely on the result of the estimation but the optimization of the correlation between the independent vector and the reference vectors. Besides, unlike standard motion correction methods, McIVA is a non-linear approach which could correct the slice and even voxel-based motion artifacts. What is more, as McIVA is applied to every volume including both large and small estimated motion, it can improve the motion correction performance from a large scale.  In the next part, we will demonstrate that, combined with nuisance regression, McIVA could achieve excellent performance on motion denoising.  5 10 15 20 25 30 35 40 45ROI Distance (voxel)00.10.20.30.40.5ROI CorrelationOriginalMCIVASPMDof6Dof1252  4.7 Comparison with common denoising methods The experiment in this part was carried out on the dataset-2 which asked for a full preprocessing pipeline. The proportion of uncorrected p < 0.05 is compared to measure the efficacy of those denoising approaches.  4.7.1 Effect of preprocessing Since years ago, brain scientists have noticed that different preprocessing methods will have a significant effect on the fMRI data (Strother, 2006). Some recent studies have pointed out that both standard and non-standard data pre-processing methods may have significant effects (negative or positive) on the quality of fMRI data analysis (Aurich, 2015). However, few of them give a functional correlation-based analysis of the effects of the preprocessing methods. Most fMRI analyses implicitly assume that the preprocessing pipeline give optimal outputs, or their results are relatively insensitive to the chosen set of preprocessing pipelines (Churchill, 2012), which is not always the case. Although we have tried to decrease the preprocessing procedures in order to relieve their impacts, some steps are still necessary for post analysis of fMRI signal. In addition to some essential procedures like slice-timing correction and realignment, the processed volumes need to be normalized to the same shared space for group analysis. The MNI152 was adopted in our experiments. In order to draw a reliable conclusion, the estimation of the effects of the preprocessing methods on the raw fMRI data was proceeded at first. As motion's effects are highly variable, motion-related signal changes are not stereotyped and dataset-wised but shared across voxels. An analysis of the individual subject was developed to get an understanding of the variety of signal changes that preprocessing methods could produce.  53   (a) Non-McIVA processed  (b) McIVA processed Figure 4-10. The effect of the preprocessing pipeline. The effect of the same preprocessing pipeline on McIVA processed and non-McIVA processed data. McIVA processed fMRI data is less sensitive to the preprocessing procedures. 54   Figure 4-10 shows the QC measures and preprocessing-caused activation on the subject level within the dataset using the toolbox CONN. These subjects were processed with the preprocessing pipeline as in Figure 3-1. The distribution of voxel-to-voxel functional correlation throughout the whole 240 scans are demonstrated in the upper figures, and the fraction of variances explained by the preprocessing is shown in the lower images. From the functional connectivity, we could figure out that when head movement occurred, two voxels may experience different kinds of disruptions. Typically, as shown in Figure 4-10 (a) and (b), the effect of motion is similar and widespread over much of the brain, leading to a positive mean value of the functional connectivity. What is more, the motion will make a more significant impact on the functional connectivity (FC) between nearby voxels under the fact that nearby voxels always experience similar disruptions, especially in resting-state fMRI when the global signal is relatively small compared with task-based fMRI. The influence of motion on the FC distribution of more distant voxels may vary from subject to subject, which could be a similar or dissimilar disruption.  It is evident that preprocessing will have a significant effect on the distribution of functional connectivity. The original data is an approximately normal distribution with a positive mean value. The preprocessing procedures will have an effect of contraction from left to right on the fMRI data, leading to a more centralized distribution which is the yellow one in Figure 4-10. The transfer of the distribution implies that false functional connectivity information is added to the fMRI data. A more intuitive result could be seen in the heat map of Figure 4-10 (a), which shows the variance explained by the preprocessing procedures of the middle slice.  55  As shown in Figure 4-10(b), McIVA-processed data appears to have a higher preprocessing resistance, which is revealed in that preprocessing procedures barely affect the distribution of the functional connectivity. Intuitively, much fewer variances were generated after the preprocessing pipeline as shown in the heat map, which resulted in the less negative or false effect on the functional connectivity. McIVA could, to a large extent, maintain the distribution of functional connectivity.  4.7.2 Comparison with denoising approaches The most noticeable advantage of McIVA was not confined to the preprocessing. In fact, combined with primary nuisances regression, McIVA could achieve much better performance on motion denoising. Two novel denoising methods aCompCor and ICA-AROMA were compared with McIVA in this part. For aCompCor, five principal components of WM and CSF were applied. For aCompCor and McIVA, they have a similar original group-wise input distribution as they used a similar preprocessing pipeline. While as ICA-AROMA asks for a specific preprocessing order, its initial distribution may be a bit different. As an ICA-based method, ICA-AROMA could automatically detect and remove the independent components of noise. As its function may overlap with CSF regression, only the realignment parameters regression is applied with ICA-AROMA here.  Figure 4-11 shows that all three methods tend to decrease the voxel-wise functional correlation. Among the three methods, ICA-AROMA appears to be less effective than the other two. Although the processed data is more centralized, the mean value of the functional connectivity is around 0.07, indicating incomplete removal of motion artifacts. Actually, from the waterfall plots of the 56  three subjects as shown in Figure 4-11(c), some motion artifacts are also visible to the naked eye after processed with ICA-AROMA. McIVA and aCompCor appear to have a comparable result on reducing the mean value of the functional connectivity which is around 0.03.  57   (a) aCompCor + 12 Realignment Parameters  (b) ICA-AROMA + 12 Realignment Parameters  (c) McIVA + CSF + WM (5p) Figure 4-11. Functional connectivity distribution of aCompCor, McIVA and ICA-AROMA. 58  a) apply McIVA for motion correction and CSF, WM regression in denoising. b) apply aCompCor (5p CSF, WM) and 12 RPs regression in denoising. c) apply ICA-AROMA and 12 RPs regression in denoising. However, apparently, McIVA has a more centralized distribution than that of aCompCor. What is more, the signal of white matter is preserved in McIVA.  In general, McIVA based method could achieve a better denoising performance than aCompCor(5P) and ICA-AROMA on the basis of preserving tissue signals. The group-based correlation analysis on the right also indicates that our approach is effective on a larger scale.  4.7.3 ROI-based functional correlation Figure 4-12 shows the functional connectivity between ROIs after motion denoising procedures with ICA-AROMA, aCompCor, and McIVA separately. For a clear demonstration, the results of 15 HC subjects are shown in the group-wise ROI-ROI correlation. The red and blue spheres represent the correlation between the 'seed' ROI and the other ROIs. Red means positive correlation and blue means a negative correlation. The lower row indicates the seed-ROI correlation of all the subjects in the dataset. Each vertical bar is proportional to the strength of correlation. Blue vertical lines represent a negative correlation.   Among the three, aCompCor appears to be less effective than the other two methods on denoising. All of the correlation values between ROIs of all the subjects are still quite high after applying aCompCor. As aCompCor only extract principal components from WM and CSF, it is less helpful on removing the motion artifacts of the other tissues. ICA-AROMA performs better than aCompCor in some cases as, not like aCompCor, ICA-AROMA is applied at the global scope. While it is not effective in all the subjects, though ICA-AROMA could be used for a single subject, 59  its effectiveness relies too much on the accuracy of the classification of the independent components.  Figure 4-12. Functional connectivity between the seed ROI and other ROIs.  The red and blue spheres represent the correlation between the 'seed' ROI and the other ROIs. Red means positive correlation and blue means a negative correlation. The lower row indicates the seed-ROI correlation of all the subjects in the dataset. Each vertical bar is proportional to the strength of correlation. Blue vertical lines represent a negative correlation.  Meanwhile, the accuracy of the classifier is limited based on the data of a specific subject. McIVA shows a strong ability to weakening the positive correlation effect caused by head motion on our dataset. Some ROIs also have a negative correlation (the blue ones), which is consistent with the ideal situation that the functional correlation should be a symmetric distribution around the lateral axis when no motion happens.  60  Chapter 5: Discussion  We have proposed a novel independent vector analysis-based motion correction method for fMRI time series which appears to have some desirable properties. Superficially, as a motion correction method, McIVA could achieve non-linear correction for motion artifacts. The implicit criterion which lies in IVA makes it possible to find the independent component in the sub-volumes which is the most related to the reference volumes. In fact, the most advantage of McIVA is to reduce the differences among the time series under the criterion of maintaining the activation information of itself.  The main benefit of McIVA is that, though based on the estimation of realignment parameters, the accuracy of McIVA is independent of the accuracy of the estimation for the reason that the estimation is only used to generate the sub-volumes rather than to correct the motion-contaminated volume directly. A big premise of all the preprocessing methods, including motion correction and denoising methods like ICA-AROMA, aCompCor as well as McIVA, is that, though contaminated by motion artifacts and some other noise, the raw fMRI time series have all the information needed by the analysis. Otherwise, we are definitely not able to recover the useful information from the raw fMRI time series. A dense search for the sub-volumes in the limited range of the possible transformation was performed first in order to obtain all the information of a specific volume. This idea was inspired by another motion correction method called Motion Corrected Independent Component Analysis (MCICA) (Liao, 2005). However, in that method, the author assumed that the actual volume was a linear combination of the sub-volumes and the optimal coefficients could be gotten from optimizing a mutual information-based cost function, which is not always the truth. 61  The actual volume hides in the space created by the sub-volumes, and most of the time the relationship between the actual volume and the sub-volumes is non-linear. That is why we made use of the properties of sorting and correlation of independent vector analysis to extract the motion-corrected volume from the sub-volume space.  What some other conventional methods like ICA-AROMA and aCompCor have in common is that they try to extract useless signals from the raw fMRI time series. However, the difficulty of fMRI analysis lies in the fact that it is hard to obtain the ground truth of the volume. That is to say, the more the components are removed, the more useful information may be eliminated as well. The criterion of McIVA is to try to include as more information as possible under the premise of maximizing the correlation among the volumes. The results in the last chapter have indicated that McIVA could achieve a better performance than the other two methods on removing the confounding effects of head motion on our dataset though there are still some remaining artifacts.  5.1 Difference between McIVA and MCICA Though MCIVA and MCIVA have some thoughts in common, the core algorithms are entirely different. MCICA argue that movement during fMRI data acquisition results in a simultaneous increase in the joint entropy of the observed time-series and a decrease in the joint entropy of a nonlinear function of the derived spatially independent components calculated by ICA. The entropy difference was proposed as a criterion for motion correction and maximized. The registered data matrix  𝑋° = ±𝑅J𝑝, 𝑥4, 𝑥; …𝑥²³ (5.1) 62   is used to maximize the entropy difference E{log |𝐽|}. Here, 𝑥¸	(𝑖 = 2,… 𝑞) represents the subset of volumes that act as the “reference volumes” for the fMRI time-series, and the motion-corrupted time point is now replaced by 𝑹𝑻𝒑, a linear sum of the basis images 𝑹 in weighted by the motion-compensation vector 𝒑.  Actually, it is a sufficient and unnecessary condition for motion correction. Maximizing the entropy difference could alleviate the interference introduced by head motion to some extent, but it does not indicate that the head motion is removed.  Besides, MCICA used 𝑋° above, assuming the base images 𝑹 and the reference volumes are in the same observation dataset and the motion-corrupted volume and the references have the same sources, ignoring the independent features of the motion-corrupted volume. The result of MCICA might be misled by a false ICA result.  While what McIVA is trying to achieve is to extract the right independent component from the subspace of the basis images 𝑹. The structure of IVA_Core in Figure 2-1 is shown in Figure 5-1.  Given the observation datasets 𝑥¸ 𝑥¸ = 	~𝑎¸¿ ∘ 𝑠¿Á¿ (5.2) the source component vector 𝑠¸ is  63  𝑠¸ ≈ 	?̂?¸ = 	~𝑤¸¿T¿ ∘ 𝑥¿ (5.3) McIVA makes use of the dependent information between the source component vectors and higher-order-statistics, and the combination of multiple datasets makes it has a better performance than conventional ICA, MCCA, and group-ICA during the decomposition. McIVA also considers that the motion-corrupted volume and the reference volumes have different but highly dependent sources which are closer to the practical situation.   Figure 5-1. Structure of the IVA Core. The IVA core is the IVA structure in red dash line in Figure 2-1.   Then we also have some experiments to compare the performance of MCICA and McIVA to demonstrate the improvement of the proposed method. First, the motion estimation of the corrected data of the two methods is shown in Figure 5-2. The results of the estimation are separated into six parameters. Compared to MCICA, McIVA shows a lower and more stable motion estimation results, especially in the rotation aspect. The estimated motion was controlled under around 64  −0.3	~	0.2	𝑚𝑚  in translation and −0.1	~	0.03	𝑑𝑒𝑔𝑟𝑒𝑒  in rotation with a voxel size of 3 × 3 × 3	𝑚𝑚, which is much lower than that of MCICA.  Figure 5-2. Motion estimation of the fMRI data motion corrected with MCICA and McIVA.  Secondly, as the cost function of MCICA is based on the entropy change H(y) – H(x), we also analyzed its robustness on real fMRI data.  0 50 100 150 200 250Translation-0.1-0.0500.050.1mmMCICAx translationy translationz translation0 50 100 150 200 250Rotation-0.2-0.100.10.2degreepitchrollyaw0 50 100 150 200 250Translation-0.1-0.0500.050.1mmMcIVAx translationy translationz translation0 50 100 150 200 250Rotation-0.2-0.100.10.2degreepitchrollyaw65   Figure 5-3. Relationship of H(y) H(x) and the estimated motion in real fMRI data. As shown in Figure 5-3, when timepoints increases and for real fMRI data it is not as clean as simulation, the relationship between the estimated motion and the entropy H(y) and H(x) is quite weak. Besides, the complexity of real fMRI data will make the cost function fall into the local optimum easily or cannot get the global optimum. 0 50 100 150 200 250Translation-0.3-0.2-0.100.10.2mmDisplacementx translationy translationz translation0 50 100 150 200 250Rotation-0.3-0.2-0.100.10.20.30.4degreepitchrollyaw0 50 100 150 200 250time points1.231.2351.241.2451.25H(x)MCICA0 50 100 150 200 250time points7.867.887.97.927.947.967.988H(y)66   (a) Cost function vs. training steps  (b) H(y) vs train steps     (c) H(x) vs train steps Figure 5-4. Convergence properties of the cost function of MCICA. The cost function change of MCICA is shown in Figure 5-4. The value converges after about 40 iterations with some fluctuations. While, actually, the convergence pattern is mainly due to the algorithm of ICA which takes the entropy H(y) as the cost function (Figure 5-4(b)) to minimize 0 10 20 30 40 50Steps6.66.626.646.666.686.76.726.746.76cost functionMCICA0 10 20 30 40 50Steps7.847.867.887.97.927.947.967.988H(y)MCICA0 10 20 30 40 50Steps1.2251.231.2351.241.2451.25H(x)MCICA67  the mutual information between the independent components but not the motion correction. Figure 5-4(c) shows that H(x) changes little with iteration. In fact, the process of convergence is also quite unstable even though the learning rate is already set to be quite small.  Figure 5-5. Convergence properties of the cost function of McIVA.  Compared with MCICA, the cost function of McIVA, as shown in Figure 5-5, tends to be quite stable under different motion conditions. Meanwhile, for the algorithm of IVA, global optimum could be validated through the permutation results.  According to references (Power) (Whitfield-Gabrieli, 2012) (Power, 2012), voxel-wise functional connectivity based on preprocessed fMRI data (where pre-processing includes standard processing, motion correction, brain segmentation and further denoising as the final step) is a common strategy to assess the data cleaning/interference removal effects of pre-processing. With 68  motion correction component being the only different component in the pre-processing pipeline, we can compare the efficiency of MCICA and McIVA. We can see that the performance of motion correction has a significant effect on the final result of denoising. Intuitively, a better motion correction contributes to a more accurate brain segmentation, and therefore may lead to a better final denoising result as indicated by a desired function connectivity distribution. Figure 5-6 is an example of voxel-wised functional connectivity of fMRI data motion corrected with MCICA and McIVA.     (a) MCICA corrected     (b) McIVA corrected Figure 5-6. Voxel-wised functional connectivity distribution on subject level.  In the fMRI literature, the motion of fMRI data could be separated as low motion (less than 1mm), high motion (1 to 3mm), severe motion (larger than 3mm) and currently no method can satisfactorily correct such large motion. In practice, subjects with this kind of large motion 69  generally should be excluded for further data analysis. Here we investigate the functional connectivity which covers all the motion conditions to compare the performance of MCICA and McIVA. Figure 5-7 shows the functional connectivity distribution of a batch of fMRI data. The distributions are a bit similar but actually quite different. The statistical data were calculated in Table 5-1. Voxel-wised functional connectivity based on different motion conditions.   Table 5-1. Voxel-wised functional connectivity based on different motion conditions. Subject Index 1 2 3 4 5 6 7 Mean Translation (mm) 0.2166 0.2552 0.2703 0.2715 0.8694 1.116 1.7328 Mean Rotation (degree) 0.2663 0.7122 0.3716 0.6411 0.5965 1.6335 1.3231 Max Translation 0.4631 0.4947 0.7136 0.4838 1.6283 1.928 4.4959 Max Rotation 0.4903 1.0446 0.8354 1.0416 1.423 2.8221 3.3081    Functional Connectivity (mean value) McIVA 0.02 0.02 0.04 0.03 0.05 0.02 0.05 MCICA 0.04 0.09 0.08 0.07 0.06 0.03 0.05   Degree of Freedom McIVA 70.6 70.6 70.6 70.6 70.6 70.6 70.6 MCICA 67.1 67.1 67.1 67.1 67.1 67.1 67.1  As shown in Table 5-1, McIVA have a better performance than MCICA on improving the functional connectivity on all the motion condition whose translation is ranging from 0 to 4.4959mm and rotation is ranging from 0 to 3.3081 degree. Only in the severe motion condition (red bold number), MCICA has a comparable value with that of McIVA but McIVA preserves more degree of freedom. Here the motion was defined as translation and rotation separately and 70  subject 7 in Table 5-1 has the largest motion in our dataset. For resting-state fMRI data, it is a considerable head motion.   (a) FC of fMRI data motion corrected by MCICA  (b) FC of fMRI data motion corrected by McIVA Figure 5-7. Distribution of functional connectivity of fMRI data on group level.  Then the distribution of subject-level QC measures of the two approached were calculated, as shown in Figure 5-8 and  Figure 5-9. We could find that, except subject 4 (the subject index here is not the same as Table 5-1), McIVA is quite stable under different kind of conditions. McIVA could preserve more valid scans of the subjects. The motion and global signal changes of McIVA 71  also have a uniform distribution. By comparison, MCICA is more sensitive to different subject conditions.   Figure 5-8. Distribution of subject-level QC measures (a).  Distribution of subject-level QC measures of the fMRI data motion corrected with McIVA. 72    Figure 5-9. Distribution of subject-level QC measures (b). Distribution of subject-level QC measures of the fMRI data motion corrected with MCICA.  Then we dig further to the subject 4 to see why s4 was detected as an outlier in both MCICA and McIVA. As shown in Figure 5-8 and Figure 5-9, tens of scans of subject 4 were classified as invalid scans after the preprocessing in both MCICA and McIVA. The estimated motion of subject 4 is shown in Figure 5-10. The range of the motion is not quite high which is from 0.5 to 1.5mm in translation and less than 1.3 degree in rotation, but actually there was a violent tremor throughout 73  the whole fMRI acquisition process which is fatal to fMRI data. In this condition that nearly all the scans are corrupted, reacquisition should be the best choice.  Figure 5-10. Motion estimation of the original fMRI data of the subject 4.  In conclusion, through several groups of comparison test, McIVA shows a better performance than MCICA under different motion and subject conditions.  5.2 Concern about white matter Most effects motion produces across the majority of voxels in the brain are non-uniform, which suggests that the variance should be effectively reduced through regressing confounds derived from white matter, gray matter, CSF or other brain voxels. The effectiveness of the method aCompCor shown in Figure 4-11 a) has demonstrated that. 74   The more confounds are regressed out, the ‘cleaner’ fMRI signal will be gotten. And therefore, an ideal correlation map will be derived. Motion parameters have been commonly used as regressors when applying GLM, and the regression of CSF is acceptable as there is no BOLD signal in CSF tissue. However, it is a controversial issue whether the WM signal should be regressed for the reason that white matter also plays a vital role in the analysis of brain connectivity. Besides, more importantly, there is no ground truth for fMRI studies, which means that every confounds regressed out could have relevant information about the brain activities.   Some recent studies have begun to analysis the bold oxygenation level dependent (BOLD) signals in white matter and attempt to detect those resting-state connectivity (Wu, 2016).   Figure 5-11. Motion denoising using McIVA combined with nuisance regression. More nuisance regressors such as WM, CSF, age and gender were regressed out in the nuisance regression.  75  Unlike voxel to voxel functional correlation, which is a global estimation and evaluation of the motion artifacts, ROI based connectivity analysis is quite relevant to the choice of specific ROIs. What we concerned is that, whether the WM should be considered as a nuisance regressor. The result of McIVA combined with regression of RP, CSF and WM are shown in Figure 5-11 and Figure 5-12.  For the voxel-to-voxel functional connectivity, a much more centralized and clean data could be gotten on the group level as shown in Figure 5-11, indicating almost complete removal of motion artifacts. While, by contrast, the result of ROI based functional connectivity on subject level as shown in Figure 5-12, has no improvement compared to the result in Figure 4-12(c), which means that for the study of some ROI based functional connectivity, regressing WM is just a cosmetic procedure and has no improvement in removing the motion artifacts.  76   Figure 5-12. Group-wise ROI-to-ROI correlation of Figure 5-11. CSF and WM were regressed out during the denoising procedure. Red means positive correlated and blue means negative correlated. The height of the bar represents the degree of the correlation.  5.3 Number of reference volumes According to the theory and the main structure of McIVA, the number of reference volumes is supposed to have a big effect on the performance of McIVA. In order to analyze this factor, we also had some experiments under a different number of reference volumes. 77   (a) 50 reference volumes   (b) 80 reference volumes  (c) 120 reference volumes Figure 5-13. ROI based correlation under different number of reference volumes of McIVA. a) use 50 as reference volumes; b) use 80 as reference volumes; c) use 120 as reference volumes.  Figure 5-13 shows the ROI-based functional correlation under a different number of reference volumes. References in McIVA should be those volumes with relatively less head motion. Although we have applied motion correction among the reference volumes first, the original quality of the reference volumes matters as well. If the number of the reference volumes is too small, the result of PCA for generating the ‘observation’ dataset will be quite limited. On the contrary, if the number is too big, like 80 or 120 shown in Figure 5-13, some volumes with large 5 10 15 20 25 30 35 40 45ROI Distance (voxel)00.050.10.150.20.250.30.350.40.450.5ROI CorrelationOriginalMCIVASPMDof6Dof125 10 15 20 25 30 35 40 45ROI Distance (voxel)00.050.10.150.20.250.30.350.40.450.5ROI CorrelationOriginalMCIVASPMDof6Dof125 10 15 20 25 30 35 40 45ROI Distance (voxel)-0.100.10.20.30.40.5ROI CorrelationOriginalMCIVASPMDof6Dof1278  head motion will be calculated as reference volumes and affect the performance of motion correction.  5.4 Drawbacks Though McIVA shows more advantages than other popular methods on our dataset, there are some drawbacks which should be improved as well. First, as the extraction of sub-volumes and operations like PCA and IVA ask for a large number of computing resources, improved optimization based on GPU and distributed system should be deployed to increase the speed of this approach. Second, the extraction of sub-volumes is isotropic. A more reasonable sampling method should be developed that more sub-volumes should be sampled in the direction of head movement.    79  Chapter 6: Conclusion and Future work  6.1 Conclusion and Contribution In this thesis, we proposed the McIVA, a novel strategy for removing motion artifacts from fMRI data. The McIVA exploits independent vector analysis to recover the brain volumes from head motion in a data-driven approach. When head motion occurs, it will not be appropriate to apply any post-processing method before sufficiently eliminating the influence of such motion. The proposed McIVA can avoid removing volumes from the raw fMRI data, and largely preserve temporal degrees of freedom. Besides, by taking advantage of the essential information of a specific subject itself, McIVA could be robustly deployed across datasets without requiring any prior heuristic knowledge, manual annotation or training. The proposed McIVA is a robust and accurate method for pre-processing of fMRI data and could lead to a more accurate analysis of functional connectivity between ROIs.  In chapter 2, we developed the framework of motion correction with independent vector analysis to solve the problem of inherent non-linearity of head motion in fMRI signals. The desired properties IVA (e.g., making use of the correlations between source component vectors) make it possible for us to extract the most correlated independent volume from the ‘observation’ volumes.  In chapter 3, we introduced the pre-processing pipeline for the proposed McIVA approach. As a motion correction method, McIVA is suggested to be applied at the beginning of the pre-processing pipeline. Realignment and reslicing procedures can be replaced by McIVA.  80  In chapter 4, we evaluated the motion correction performance of the proposed McIVA and compared its improvement on denoising with two popular methods, aCompCor and ICA-AROMA. The proposed McIVA shows improved results.   In chapter 5, we discussed several factors that have an effect on the results of pre-processing in related fMRI analysis. As the proposed McIVA was inspired by MCICA, the differences between the two methods were investigated from both the theoretical implementation and experiments. Besides, we also investigated the influences of the number of reference volumes and discussed some current drawbacks in McIVA.   From our study, the proposed McIVA is suggested to replace the conventional motion correction and realignment steps in the fMRI pre-processing pipeline. Some recent studies point out that the interpolation operation in reslicing and realignment will have a significant effect on motion estimation (Power, 2017). In fact, from the aspect of fMRI analysis, interpolation is employed mainly for its computational convenience, while its meaningfulness in practice may not be true. When applying the interpolation operation, the intensity values of two adjacent voxels are assumed to be continuous which is actually not always the case.  By contrast, the proposed approach is more practically meaningful when applied for motion correction. The results in chapter 4 have demonstrated that the proposed McIVA could achieve a lower MAE than that of the interpolation. Meanwhile, when combined with McIVA, the confounds of CSF are also recommended to be regressed out in the denoising procedure (Geissler, 2005). However, for the confounds in the white matter, we recommend that whether regressing out those confounds or not should depend on what kind of analysis to be taken in the post-processing for different purposes. 81   6.2 Future Considerations In Chapter 5, we have discussed the results and some current drawbacks in applying the proposed McIVA. Focusing on overcoming these shortages, we plan to propose a few improvements on the current McIVA approach to provide a more convenient and reliable tool for motion correction of fMRI data.  Parallel computing will be employed to accelerate the process of preprocessing in McIVA. The accuracy of McIVA is based on a large number of computations including generating sub-volumes, utilizing PCA for generating the ‘observation’ dataset of the motion-corrupted volume and applying IVA to extract the corresponding component for motion correction. GPU acceleration is also a solution for improving the matrix related operations in McIVA.  Besides, the choice of the number of reference volumes in McIVA is still quite heuristic. McIVA shows excellent performances on our PD and HC datasets; while when applying it to other datasets, the number of reference volumes may be different. We recommend setting the number of the reference volumes to be about 15% to 20% of the total time points. Due to the lack of experimental data currently available, we did not explore further in this aspect. We plan to further study this and make the choice automatic in the future.  Furthermore, we are working to compile the proposed algorithm into a toolbox or a library for its integration into other fMRI preprocessing tools.  82  Bibliography  Aurich, N. K. e. a., 2015. Evaluating the reliability of different preprocessing steps to estimate graph theoretical measures in resting state fMRI data. Frontiers in neuroscience, p. 48. Behzadi, Y. e. a., 2015. A component based noise correction method (CompCor) for BOLD and perfusion based fMRI. Neuroimage, pp. 90-101. Caballero-Gaudes, C. a. R. C. R., 2017. Methods for cleaning the BOLD fMRI signal. Neuroimage, pp. 128-149. Churchill, N. W. e. a., 2012. Optimizing preprocessing and analysis pipelines for single‐subject fMRI. I. Standard temporal motion and physiological noise correction methods. Human brain mapping, pp. 609-627. Crinion, J. e. a., 2007. Spatial normalization of lesioned brains: performance evaluation and impact on fMRI analyses. Neuroimage, Volume 37.3, pp. 866-875. Friston, K. J. e. a., 1996. Movement‐related effects in fMRI time‐series. Magnetic resonance in medicine, pp. 346-355. Gawryluk, J. R. E. L. M. a. R. C. D., 2014. Does functional MRI detect activation in white matter? A review of emerging evidence, issues, and future directions. Frontiers in neuroscience, p. 239. Geissler, A. e. a., 2005. Influence of fMRI smoothing procedures on replicability of fine scale motor localization. Neuroimage, Volume 24.2, pp. 323-331. Grootoonk, S. e. a., 2000. Characterization and correction of interpolation effects in the realignment of fMRI time series. NeuroImage, Volume 11.1, pp. 49-57. 83  Hutton, C. e. a., 2011. The impact of physiological noise correction on fMRI at 7 T. Neuroimage, pp. 101-112. Jenkinson, M. e. a., 2002. Improved optimization for the robust and accurate linear registration and motion correction of brain images. Neuroimage, Volume 17.2, pp. 825-841. Kim, T. T. E. a. T.-W. L., 2006. Independent vector analysis: An extension of ICA to multivariate components. Berlin, Heidelberg, Springer. Liao, R. J. L. K. a. M. J. M., 2005. An information-theoretic criterion for intrasubject alignment of FMRI time series: motion corrected independent component analysis. IEEE Transactions on Medical Imaging, pp. 29-44. Lund, T. E. e. a., 2005. Motion or activity: their role in intra-and inter-subject variation in fMRI. Neuroimage, pp. 960-964. M.W. Woolrich, S. J. B. P. M. C. S. M. T. B. C. B. M. J. S. S., 2009. Bayesian analysis of neuroimaging data in FSL. NeuroImage. Penny, W. D. e. a., 2011. eds. Statistical parametric mapping: the analysis of functional brain images. Elsevier. Poldrack, R. A., 2007. Region of interest analysis for fMRI. Social cognitive and affective neuroscience, pp. 67-70. Power, J. D. B. L. S. a. S. E. P., 2015. Recent progress and outstanding issues in motion correction in resting state fMRI. Neuroimage, pp. 536-551. Power, J. D. e. a., 2012. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage, pp. 2142-2154. 84  Power, J. D. e. a., 2012. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage, 59(3), pp. 2142-2154. Power, J. D. e. a., 2013. Steps toward optimizing motion artifact removal in functional connectivity MRI; a reply to Carp. Neuroimage. Power, J. D. e. a., 2014. Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage, pp. 320-341. Power, J. D. e. a., 2014. Methods to detect, characterize, and remove motion artifact in resting state fMRI. Neuroimage, pp. 320-341. Power, J. D. e. a., 2017. Temporal interpolation alters motion in fMRI scans: Magnitudes and consequences for artifact detection. PloS one, Volume 12.9. Pruim, R. H. e. a., 2015. Evaluation of ICA-AROMA and alternative strategies for motion artifact removal in resting state fMRI. Neuroimage, pp. 278-287. Pruim, R. H. e. a., 2015. ICA-AROMA: A robust ICA-based strategy for removing motion artifacts from fMRI data. Neuroimage, pp. 267-277. Satterthwaite, T. D. e. a., 2013. An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage, pp. 240-256. Satterthwaite, T. D. e. a., 2013. An improved framework for confound regression and filtering for control of motion artifact in the preprocessing of resting-state functional connectivity data. Neuroimage, pp. 240-256. 85  Strother, S. C., 2006. Evaluating fMRI preprocessing pipelines. IEEE Engineering in Medicine and Biology Magazine, pp. 27-41. Tsai, A. e. a., 1999. Analysis of functional MRI data using mutual information. Berlin, Heidelberg,, Springer. Van Dijk, K. R. M. R. S. a. R. L. B., 2012. The influence of head motion on intrinsic functional connectivity MRI. Neuroimage, Volume 59.1, pp. 431-438. Whitfield-Gabrieli, S. a. A. N.-C., 2012. Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks. Brain connectivity, pp. 125-141. Whitfield-Gabrieli, S. a. A. N.-C., 2012. Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks. Brain connectivity, 2(3), pp. 125-141. WhitfieldWhitfield-Gabrieli, S. a. A. N.-C., 2012. Conn: a functional connectivity toolbox for correlated and anticorrelated brain networks. Brain connectivity, pp. 125-141. Wu, T.-L. e. a., 2016. Effects of anesthesia on resting state BOLD signals in white matter of non-human primates. Magnetic resonance imaging, pp. 1235-1241. Yan, C.-G. e. a., 2013. A comprehensive assessment of regional variation in the impact of head micromovements on functional connectomics. Neuroimage, pp. 183-201.   

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0379472/manifest

Comment

Related Items