Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Multimodal fusion for assessing functional segregation and integration in the human brain Yoldemir, Burak 2016

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

Item Metadata


24-ubc_2016_may_yoldemir_burak.pdf [ 5.17MB ]
JSON: 24-1.0225963.json
JSON-LD: 24-1.0225963-ld.json
RDF/XML (Pretty): 24-1.0225963-rdf.xml
RDF/JSON: 24-1.0225963-rdf.json
Turtle: 24-1.0225963-turtle.txt
N-Triples: 24-1.0225963-rdf-ntriples.txt
Original Record: 24-1.0225963-source.json
Full Text

Full Text

Multimodal Fusion for AssessingFunctional Segregation andIntegration in the Human BrainbyBurak YoldemirB.Sc., Bogazici University, 2009M.Sc., Bogazici University, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)February 2016c© Burak Yoldemir 2016AbstractMental and neurological diseases account for a major portion of the globaldisease burden. Neuroimaging has greatly contributed to the characteriza-tion and understanding of such disorders by enabling noninvasive investi-gation of the human brain. Among neuroimaging technologies, magneticresonance imaging (MRI) stands out as a relatively widespread and safeimaging modality that can be sensitized to capture different aspects of thebrain. Historically, MRI studies have investigated anatomy or function ofthe brain in isolation, which created an apparent dichotomy. In this thesis,we aim to bridge this divide using novel multimodal techniques. In partic-ular, we present techniques to reconcile information regarding anatomicaland functional connectivity (AC and FC) in the brain estimated from dif-fusion MRI (dMRI) and functional MRI (fMRI) data, respectively. Ourfirst contribution is to show that the consistency between AC and FC isunderstated when standard analysis methods are used. We illustrate howthe estimation of AC can be improved to increase the AC-FC consistency,which facilitates a more meaningful fusion of these two types of information.Specifically, we propose to improve AC estimation by the use of a dictio-nary based super-resolution approach to increase the spatial resolution indMRI, reconstructing the white matter tracts using global tractography in-stead of conventional streamline tractography, and quantifying AC usingfiber count as the metric. Our second contribution is to develop novel mul-timodal approaches for investigating functional segregation and integrationin the human brain. We show that task fMRI data can be fused with dMRIand resting state fMRI data to mitigate the effects of noise and deconfoundthe effects of spontaneous fluctuations in brain activity on activation detec-tion. Further, we show that sensitivity in unraveling the modular structureof the brain can be increased by fusing dMRI and fMRI data. Our resultscollectively suggest that combining dMRI and fMRI data outperforms clas-sical unimodal analyses in understanding the brain’s organization, bringingus one step closer to understanding the most complex organ in the humanbody.iiPrefaceThe research presented herein was approved by the UBC Clinical ResearchEthics Board, certificate number: H15-02343. This thesis is primarily basedon the following papers, resulting from collaboration between multiple re-searchers.Studies described in Chapter 2 have been published in:[P1] B. Yoldemir, B. Ng, R. Abugharbieh. Effects of tractography ap-proach on consistency between anatomical and functional connectiv-ity estimates. In Proceedings of the IEEE International Symposiumon Biomedical Imaging (ISBI’14), pp. 250-253, Beijing, China, April2014. Oral presentation - Acceptance rate: ∼20%.[P2] B. Yoldemir, M. Bajammal, R. Abugharbieh. Dictionary based super-resolution for diffusion MRI. In Proceedings of the MICCAI Workshopon Computational Diffusion MRI, pp. 194-204, Boston, MA, Septem-ber 2014.[P3] M. Bajammal, B. Yoldemir, R. Abugharbieh. Comparison of struc-tural connectivity metrics for multimodal brain image analysis. In Pro-ceedings of the IEEE International Symposium on Biomedical Imaging(ISBI’15), pp. 934-937, Brooklyn, NY, USA, April 2015.Studies described in Chapter 3 have been published in:[P4] B. Yoldemir, B. Ng, R. Abugharbieh. Deconfounding the effects ofresting state activity on task activation detection in fMRI. In Proceed-ings of the MICCAI Workshop on Multimodal Brain Image Analysis,LNCS, vol. 7509, pp. 51-60, Nice, France, October 2012.[P5] B. Yoldemir, B. Ng, T. Woodward, R. Abugharbieh. Fiber connectiv-ity integrated brain activation detection. In Proceedings of the 23rdBiennial International Conference on Information Processing in Medi-cal Imaging (IPMI’13), LNCS, vol. 7917, pp. 135-146, Asilomar, CA,June 2013. Acceptance rate: 32.1%.iiiPrefaceStudies described in Chapter 4 have been published in:[P6] B. Yoldemir, B. Ng, R. Abugharbieh. Overlapping replicator dynam-ics for functional subnetwork identification. In Proceedings of the In-ternational Conference on Medical Image Computing and ComputerAssisted Intervention (MICCAI’13), LNCS, vol. 8150, pp. 682-689,Nagoya, Japan, September 2013. Acceptance rate: 32.8%.[P7] B. Yoldemir, B. Ng, R. Abugharbieh. Coupled stable overlappingreplicator dynamics for multimodal brain subnetwork identification.In Proceedings of the 24th Biennial International Conference on Infor-mation Processing in Medical Imaging (IPMI’15), LNCS, vol. 9123,pp. 770-781, Isle of Skye, Scotland, June 2015. Acceptance rate:32.3%.[P8] B. Yoldemir, B. Ng, R. Abugharbieh. Stable overlapping replicatordynamics for brain community detection. IEEE Transactions on Med-ical Imaging, vol. 35, no. 2, pp. 529-538, February 2016. Impactfactor: 3.39, h5-index: 63.[P1, P4, P6-8]: I was the primary author of these papers and main con-tributor to the design, implementation and testing of the proposed methodsunder the supervision of Dr. Rafeef Abugharbieh. Dr. Bernard Ng helpedwith his valuable input in improving the methodology and design of theexperiments. The papers were edited by all co-authors.[P2]: I contributed the application idea conception, preprocessing of thedata and parcellation, and the validation scheme. Mohammad Bajammalcontributed the algorithmic idea conception, implementation of the methodand validation scheme, as well as generation of the results. I contributedthe majority of the effort for manuscript writing. The paper was edited byall co-authors.[P3]: I contributed the paper idea, preprocessing of the data and parcel-lation, and implementation of the tractography method and three of thefour anatomical connectivity metrics. Mohammad Bajammal contributedthe implementation of the fourth anatomical connectivity metric and thevalidation code as well as generated the results. I and Mohammad Bajam-mal contributed equally to writing the paper. The paper was edited by allco-authors.ivPreface[P5]: I was the primary author of this paper and main contributor to thedesign, implementation and testing of the proposed method under the su-pervision of Dr. Rafeef Abugharbieh. Dr. Bernard Ng assisted in improvingthe methodology and design of the experiments. Dr. Todd Woodward pro-vided data and clinical input. The paper was edited by all co-authors.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . 31.2.1 Objectives . . . . . . . . . . . . . . . . . . . . . . . . 31.2.2 Research Questions Addressed . . . . . . . . . . . . . 41.3 Brain Imaging for Macroscale Connectomics . . . . . . . . . 51.3.1 Diffusion Magnetic Resonance Imaging . . . . . . . . 51.3.2 Functional Magnetic Resonance Imaging . . . . . . . 61.4 Anatomical and Functional Brain Connectivity . . . . . . . . 71.4.1 Estimating Anatomical Connectivity . . . . . . . . . 71.4.2 Estimating Functional Connectivity . . . . . . . . . . 121.4.3 Relationship between Anatomical and Functional BrainConnectivity . . . . . . . . . . . . . . . . . . . . . . . 131.5 Functional Segregation in the Brain: Activation Detection . 151.5.1 Unimodal Approaches . . . . . . . . . . . . . . . . . . 151.5.2 Multimodal Approaches . . . . . . . . . . . . . . . . . 16viTable of Contents1.6 Functional Integration in the Brain: Subnetwork Identifica-tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.6.1 Unimodal Approaches . . . . . . . . . . . . . . . . . . 171.6.2 Multimodal Approaches . . . . . . . . . . . . . . . . . 191.7 Thesis Contributions . . . . . . . . . . . . . . . . . . . . . . 201.7.1 Anatomical and Functional Brain Connectivity . . . . 201.7.2 Multimodal Brain Activation Detection . . . . . . . . 211.7.3 Multimodal Brain Subnetwork Identification . . . . . 212 Anatomical and Functional Brain Connectivity . . . . . . . 222.1 Super-resolution for dMRI . . . . . . . . . . . . . . . . . . . 242.1.1 Acquisition Model . . . . . . . . . . . . . . . . . . . . 252.1.2 Dictionary Construction . . . . . . . . . . . . . . . . 252.1.3 Super-Resolved Volume Generation . . . . . . . . . . 262.1.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . 272.1.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 302.2 Estimating Anatomical Connectivity . . . . . . . . . . . . . . 332.2.1 Modelling of Local Diffusion Properties . . . . . . . . 342.2.2 Tractography . . . . . . . . . . . . . . . . . . . . . . . 342.2.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . 352.2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 372.3 Quantifying Anatomical Connectivity . . . . . . . . . . . . . 382.3.1 Anatomical Connectivity Metrics . . . . . . . . . . . 392.3.2 Experiments . . . . . . . . . . . . . . . . . . . . . . . 402.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 422.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433 Multimodal Brain Activation Detection . . . . . . . . . . . 443.1 Modelling Intrinsic Functional Connectivity . . . . . . . . . . 453.1.1 Group Activation Detection . . . . . . . . . . . . . . 463.1.2 RS Subnetwork Detection . . . . . . . . . . . . . . . . 473.1.3 RS Activity Estimation and Removal . . . . . . . . . 483.1.4 Experiments . . . . . . . . . . . . . . . . . . . . . . . 503.1.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 523.2 Modelling Anatomical Connectivity . . . . . . . . . . . . . . 533.2.1 Random Walker for Activation Estimation . . . . . . 543.2.2 Anatomical Connectivity Prior Estimation . . . . . . 553.2.3 Activation Likelihood Estimation . . . . . . . . . . . 563.2.4 Group Activation Inference . . . . . . . . . . . . . . . 57viiTable of Contents3.2.5 Re-estimation of Group Activation Using Partial Con-nectome . . . . . . . . . . . . . . . . . . . . . . . . . 583.2.6 Experiments . . . . . . . . . . . . . . . . . . . . . . . 583.2.7 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 633.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634 Multimodal Brain Subnetwork Identification . . . . . . . . 644.1 Replicator Dynamics . . . . . . . . . . . . . . . . . . . . . . 654.2 Coupled Replicator Dynamics . . . . . . . . . . . . . . . . . 664.3 Coupled Stable Replicator Dynamics . . . . . . . . . . . . . 684.3.1 Graph Incrementation . . . . . . . . . . . . . . . . . . 684.3.2 Stability Selection . . . . . . . . . . . . . . . . . . . . 684.4 Coupled Overlapping Stable Replicator Dynamics . . . . . . 694.5 Parameter Selection and Implementation . . . . . . . . . . . 704.6 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 714.6.1 Synthetic Data . . . . . . . . . . . . . . . . . . . . . . 714.6.2 Real Data . . . . . . . . . . . . . . . . . . . . . . . . 724.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 724.7.1 Synthetic Data Results . . . . . . . . . . . . . . . . . 724.7.2 Real Data Results . . . . . . . . . . . . . . . . . . . . 734.8 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 775 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81viiiList of Tables2.1 Comparison of fiber length bias for AC metrics. . . . . . . . . 42ixList of Figures1.1 Projected average number of years in full health lost due todisability and premature death of the 2010 to 2020 birth co-hort in Canada. . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Projected indirect economic costs of selected neurological con-ditions due to working-age death and disability. . . . . . . . . 31.3 Sample workflow of a dMRI experiment. . . . . . . . . . . . . 101.4 Depiction of the false negative problem in AC estimation. . . 111.5 Sample workflow of an fMRI experiment. . . . . . . . . . . . 142.1 Qualitative comparison between the tract-density maps andfiber tracts reconstructed from the original and super-resolveddMRI data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.2 AC-FC correlation for 10 subjects with AC estimated fromthe data at its original resolution (2 mm isotropic), and high-resolution data (1 mm isotropic) obtained using trilinear in-terpolation, spline interpolation, CLASR and the proposedmethod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3 Number of inter-parcel tracts reconstructed from the data atits original resolution (2 mm isotropic), and high-resolutiondata (1 mm isotropic) obtained using trilinear interpolation,spline interpolation, CLASR and the proposed method. . . . 322.4 Fiber tracts reconstructed for a representative subject usingstreamline tractography on diffusion tensors, streamline trac-tography on ODFs, and global tractography on ODFs withan ROI placed in the corpus callosum. . . . . . . . . . . . . . 372.5 Correlation between AC and FC using four commonly usedAC metrics. FC was estimated from RS fMRI data, thusdepicts intrinsic connections. . . . . . . . . . . . . . . . . . . 412.6 Correlation between AC and FC using four commonly usedAC metrics. FC was estimated from task fMRI data acquiredduring 7 different tasks. . . . . . . . . . . . . . . . . . . . . . 41xList of Figures3.1 Graphical depiction of proposed RS activity removal approach. 473.2 Dominant principal component extracted from time coursesof intrinsically-connected brain areas. . . . . . . . . . . . . . . 493.3 Comparison of activation detection results with and withoutthe proposed RS activity removal method. . . . . . . . . . . . 513.4 Detected activation patterns with and without the proposedRS activity removal method. . . . . . . . . . . . . . . . . . . 523.5 Overview of the proposed multimodal task activation detec-tion approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.6 Quantitative results of the proposed multimodal activationdetection technique. . . . . . . . . . . . . . . . . . . . . . . . 613.7 Qualitative results of the proposed multimodal activation de-tection method. . . . . . . . . . . . . . . . . . . . . . . . . . . 624.1 Subnetwork identification accuracy on synthetic data. . . . . 734.2 Subnetwork identification results on real data. . . . . . . . . . 754.3 Visuospatial subnetwork identified by CSORD at the grouplevel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.4 Right executive control subnetwork identified by CSORD atthe group level. . . . . . . . . . . . . . . . . . . . . . . . . . . 764.5 A group subnetwork identified by CSORD that comprises thevisual corticostriatal, CTC, and STC loops. . . . . . . . . . . 76xiList of AbbreviationsAC Anatomical connectivityBOLD Blood oxygen level dependentCIS Connected Iterative ScanCRD Coupled Replicator DynamicsCSORD Coupled Stable Overlapping Replicator DynamicsCSRD Coupled Stable Replicator DynamicsCSA Constant solid angleCTC Cerebello-thalamo-corticalDALY Disability-adjusted life yearsdMRI Diffusion magnetic resonance imagingDSI Diffusion spectrum imagingDTI Diffusion tensor imagingDW Diffusion weightedEM Expectation-maximizationFA Fractional anisotropyFC Functional connectivityFDR False discovery ratefMRI Functional magnetic resonance imagingFP False positiveGBD Global Burden of DiseasesGGG Gamma-Gaussian-GammaGLM General linear modelGMM Gaussian mixture modelHARDI High angular resolution diffusion imagingHCP Human Connectome ProjectHRF Hemodynamic response functionICA Independent component analysisICC Intraclass correlationLSC Left somatomotor cortexMNI Montreal Neurological InstituteMRF Markov random fieldxiiList of AbbreviationsNC Normalized cutODF Orientation distribution functionOSLOM Order Statistics Local Optimization MethodPCA Principal component analysisPCC Posterior cingulate cortexPDD Principal diffusion directionPDF Probability density functionRD Replicator DynamicsROI Region of interestRW Random walkerQBI Q-ball imagingRS Resting stateRSC Right somatomotor cortexSNR Signal-to-noise ratioSTC Striato-thalamo-corticalSVD Singular value decompositionTDI Tract-density imagingTP True positiveWHO World Health OrganizationxiiiAcknowledgementsI would like to thank my advisor, Dr. Rafeef Abugharbieh, for supervisingme throughout my doctoral studies. I owe particular thanks to Dr. BernardNg for his guidance and dedication to research. Furthermore, I would liketo thank my friends inside and outside of BiSICL for their support.Last but not least, I would like to thank my better half, Bas¸ak O¨ztas¸Yoldemir, for always being there for me.xivDedicationTo Bas¸ak, my amazing wifexvChapter 1IntroductionThe human brain is the control center of the nervous system and regulatesalmost every aspect of human behavior. It is the most complex organ inthe human body with hundreds of billions of neurons, which communicatewith each other via hundreds of trillions of connections. As such, the hu-man brain has been described as a three-pound universe [1]. Anatomically,it is an elaborate network of processing centers connected through neuralaxons. Functionally, the brain comprises specialized regions that work to-gether to enable complex cognitive processes. It is thus important to notonly characterize anatomy and function in the brain, but also interrelations,or connectivity.1.1 MotivationAccording to the Mapping Connections report published by the Governmentof Canada and Neurological Health Charities Canada in September 2014 [2],Canadians born in the current decade (2010 to 2020) who will suffer fromneurological disorders are projected to lose tens of years of life living in fullhealth. The average number of years in full health lost due to disabilityor premature death for five neurological conditions is given in Fig. 1.1. Inaddition to their physical, emotional and cognitive impacts, such disordershave immense economic costs. The Mapping Connections study suggeststhat the indirect economic cost of neurological disorders to the Canadiansociety is billions of dollars every year due to working-age death and dis-ability. Moreover, these costs are expected to steadily increase with the lifeexpectancy of the population increasing (Fig. 1.2). Worldwide impacts ofneurological conditions are becoming more widely recognised. Accordingto the latest Global Burden of Diseases (GBD) study published in Decem-ber 2012 [3], mental and neurological diseases account for a major portionof the global disease burden. To quantify the number of years lost as aresult of premature death and disability, the GBD study defines a metricdubbed disability-adjusted life years (DALYs), where one DALY equals onelost year of healthy life. In 2010, mental and neurological diseases resulted11.1. Motivation0 10 20 30 40 50Alzheimer's and other dementiasCerebral palsyEpilepsyMultiple sclerosisParkinson’s diseaseAverage number of years in full health lostWomen MenFigure 1.1: Projected average number of years in full health lost due todisability and premature death of the 2010 to 2020 birth cohort in Canada[2]. Difference between the life expectancy and the health-adjusted life ex-pectancy for each of the selected neurological conditions is 201 million DALYs with an increase of 40% since 1990 [3].During the recent years, neuroimaging has greatly contributed to thecharacterization and understanding of mental and neurological diseases.Several imaging modalities have been explored, each delineating the brainfrom a different view. Traditionally, these modalities are individually con-sidered when making inferences about the brain. The main premise of thisthesis is that fusing data from different imaging modalities should prove ben-eficial by leveraging their respective strengths. Specifically, our focus is oncombining information regarding brain activity, and the structural substratethat enables brain functions to arise.Among available imaging modalities, we chose to employ magnetic reso-nance imaging (MRI) since it is relatively widespread and safe. Global pop-ularity of MRI makes it advantageous in terms of uptake, i.e. the likelihoodof the developed techniques being adopted is high. To investigate brain func-tion, we use functional MRI (fMRI) in this thesis. fMRI has shown promisein a number of clinical applications including localization of specific brainfunctions to guide neurosurgery [4], characterization of mechanisms of dis-ease [5], prediction of treatment response [6], and characterization of disease21.2. Problem Statement0.0 0.5 1.0 1.5 2.0 2.5 3.0Alzheimer's and other dementiasCerebral palsyEpilepsyMultiple sclerosisParkinson’s diseaseTotal indirect cost ($ billions)2011 2021 2031Figure 1.2: Projected indirect economic costs of selected neurological condi-tions due to working-age (age 15 to 74 years) death and disability [2]. Dataare expressed in Canadian dollars.risk [7]. Although fMRI provides us with an effective tool to probe brain’sfunctional organization, integrating this information with anatomical con-nectivity, which essentially facilitates the functional interactions [8], mightprovide richer insights about inner dynamics of the brain. Hence in thisthesis, we propose complementing analyses of brain function with anatomi-cal connectivity information learned through diffusion MRI (dMRI), whichis currently the only modality that allows us to infer fiber pathways in thebrain in vivo and noninvasively [9].1.2 Problem Statement1.2.1 ObjectivesThe overarching objective of this thesis is to increase the sensitivity in theassessment of functional brain organization. Accurate mapping of functionalbrain organization is immensely important since it may ultimately enableidentifying reliable biomarkers for various neurological diseases. However,there are several challenges to this problem that are commonly encounteredin practice, such as poor data quality, low spatial and temporal resolu-tion, and low number of subjects. We focus on multimodal frameworks to31.2. Problem Statementmitigate such challenges. Our rationale is that combining information fromdifferent neuroimaging modalities should prove beneficial by leveraging theirrespective strengths and increasing the amount of information available.1.2.2 Research Questions AddressedAnatomical and functional connectivity (AC and FC) analyses provide twocomplementary views of brain circuitry. Motivated by how anatomicalwiring of the brain facilitates and shapes functional interactions betweenbrain regions, we propose multimodal fusion methods to integrate informa-tion regarding the two types of brain connectivity. To facilitate meaningfuldata fusion, we first address the important problem of brain connectivityestimation and quantification. Specifically, we focus on various steps of ACestimation from dMRI data and present ways of increasing the consistencybetween estimates of AC and FC. Next, we propose novel multimodal in-tegration methods to investigate the brain’s organization. It is known thatfunctional organization in the brain is governed by a subtle balance betweenfunctional segregation and integration [10]. The segregation principle statesthat some functional processes specifically engage localized and specializedbrain regions, whereas the integration principle emphasizes the importanceof connections across distributed brain regions in the emergence of brainfunction [11]. These two views of functional organization are generally ad-dressed by investigating:• which brain regions are activated upon stimulus, i.e. activation detec-tion,• which brain regions work together at rest or in response to stimulus,i.e. subnetwork (also known as community or module) identification.We propose novel methods combining fMRI and dMRI data to address thesetwo key problems.In the next sections, we first provide an overview of dMRI and fMRIand describe how they enable in vivo estimation of AC and FC. We thenprovide an outline of the two key research questions regarding functionalbrain organization, describe unimodal and multimodal techniques employedto tackle these two problems, and conclude with a summary of the thesiscontributions.41.3. Brain Imaging for Macroscale Connectomics1.3 Brain Imaging for Macroscale ConnectomicsConnectomics is the study of the human “connectome”, a term which wasintroduced in 2005 to refer to the complete map of neural connections in thebrain [12]. Recently, the term connectome is also used in a broader sense torefer to brain regions, their anatomical connections, and their functional in-teractions [13]. The connectome can be analyzed at multiple scales from themicroscale, which focuses on individual neurons and their synaptic connec-tions, to the macroscale, which focuses on large brain regions with cognitiveand behavioral associations. Currently, macroscale connectomics is by farthe most common since noninvasive tools for in vivo imaging of the humanconnectome are only available at the macroscale. Among the modalitiesused for macroscale connectomics, dMRI and fMRI are especially popular[13].1.3.1 Diffusion Magnetic Resonance ImagingdMRI is currently the only imaging modality that allows the assessmentof white matter connectivity by examining the translational displacementof water molecules [14]. Water molecules are in constant thermal motionconstrained by physical boundaries. In fibrous tissues such as axons inthe white matter, water diffuses more rapidly in the direction aligned withthe fibrous structure and more slowly perpendicular to it. Measurementsof this anisotropic diffusion hence reveal micro-structural properties of theunderlying tissue.In practice, several so-called diffusion weighted (DW) images which aresensitized to different diffusion directions are acquired, where each voxelintensity reflects the estimate of the rate of local water diffusion for thecorresponding diffusion direction. These estimates are then combined ina reconstruction step to estimate the 3D diffusion probability density func-tion (PDF) at each voxel [15]. Depending on how reconstruction is done, theimaging technique takes on different names such as diffusion tensor imaging(DTI), diffusion spectrum imaging (DSI) and Q-ball imaging (QBI) [16].Among these, DTI is the simplest and the most common imaging technique,where the local diffusion at each voxel is characterized by a 3×3 symmetric,positive semi-definite, second-order tensor called the diffusion tensor. Theinherent assumption here is that the underlying diffusion process is Gaus-sian, which is the main limiting factor of DTI. Common to all mentionedPDF reconstruction techniques, tractography is one of the most powerfultools dMRI offers. In its simplest sense, tractography is a way of delineating51.3. Brain Imaging for Macroscale Connectomicsthe fiber tracts in the brain by integrating pathways of maximum diffusioncoherence [16]. Hence, tractography allows us to infer AC patterns of thebrain in vivo and noninvasively.The most common artifacts in DW images are subject’s head motion andeddy currents arising from rapid switching of gradient pulses during dataacquisition. Eddy current artifacts are manifested in three different forms:contraction or dilation of the image, overall shift and shear [17]. Theseartifacts can be corrected with an affine registration to the b=0 image, i.e.the volume with no diffusion weighting. Head motion can be corrected witha rigid body registration to the b=0 image. These two artifacts are typicallyaccounted for using a single registration step [18].1.3.2 Functional Magnetic Resonance ImagingfMRI is an in vivo and noninvasive imaging modality for studying brainactivity, most popularly through blood oxygen level dependent (BOLD)contrast [19]. Due to hemodynamic response, i.e. rapid increase in bloodsupply to active brain regions, changes in oxygen concentration are indirectimplications of brain activity. Exploiting this mechanism, BOLD contrastimaging is based on MR images made sensitive to changes in the oxygena-tion state of hemoglobin molecule [19]. Hemoglobin behaves as a diamag-netic substance when it is saturated with oxygen (oxyhemoglobin) and asa paramagnetic substance when some oxygen molecules have been removed(deoxyhemoglobin). These changes in magnetic properties are manifestedin the fMRI observations, areas with high concentration of oxyhemoglobingiving a higher signal than areas with low concentration [20]. Broadly, thereare two different types of fMRI experiments. In the first one, the subject in-side the scanner performs a series of cognitive tasks while BOLD images arebeing collected. This experiment is typically followed by a statistical anal-ysis of the data in which signal intensity changes are compared to the taskparadigm [21]. In the second type, the subject is asked not to be engagedin any cognitive task without falling asleep, where spontaneous fluctuationsin the brain are recorded [22]. These two approaches are referred to astask fMRI and resting state fMRI (RS fMRI), and they are regarded as twodifferent modalities in this thesis.The inherently low signal-to-noise ratio (SNR) of BOLD signal (typi-cally between 0.2 and 0.5 [23, 24]) together with confounds such as scannerdrifts and subject’s head motion pose major challenges to the interpreta-tion of fMRI data. As such, fMRI data typically undergo several steps ofpreprocessing before any analysis of the data. Standard preprocessing steps61.4. Anatomical and Functional Brain Connectivityinclude slice timing correction [25], head motion correction [26], and tempo-ral autocorrelation correction [27, 28]. Slice timing correction accounts forthe fact that BOLD signals at different slices of the brain are recorded atslightly different time points. To correct for this, individual slices are tem-porally realigned to a reference slice based on their relative timing. Similarto dMRI, head motion in fMRI acquisitions is corrected by a rigid bodyregistration to a reference volume (typically the first volume). Temporalautocorrelations arise from scanner drifts, cardiac noise, respiratory noise,and head motions [28]. Since the residual of the voxel time courses is oftenassumed to be identical independently distributed (i.i.d.) in typical analysismethods, pre-whitening is performed to satisfy this assumption. In additionto these steps, spatial normalization to a common space, such as the Mon-treal Neurological Institute (MNI) space [29], can be performed to ensurespatial correspondence between different subjects [26]. There exist severaltoolboxes such as SPM [30], FSL [31] and DPARSF [32] for streamlinedpreprocessing of fMRI data.1.4 Anatomical and Functional BrainConnectivity1.4.1 Estimating Anatomical ConnectivityDiffusion Modelling and TractographyTractography using dMRI data relies on the fundamental assumption thatwater diffusion reflects the axonal orientations. The first attempt to recon-struct white matter bundles throughout the brain involved the use of dif-fusion tensors, where a 3×3 symmetric, positive semi-definite tensor field isused to model 3D Gaussian distribution of diffusing water molecules [14, 33].It was shown that the direction of least hindrance to diffusion, i.e. the prin-cipal diffusion direction (PDD), is correlated with the major eigenvector ofthe diffusion tensor at each voxel [34]. This led to several “streamline” trac-tography approaches where the main idea is to integrate the PDD at eachvoxel, using it as a surrogate for axons’ orientation [35–37]. Streamline trac-tography performed on diffusion tensors currently remains to be the normin estimating AC, mainly due to the intuitive nature of this approach: Dif-fusion tensors can be visualized as ellipsoids oriented along their PDDs, andstreamline tractography boils down to connecting the ellipsoids “pointing”at each other. However, it is now clear that PDD may not always be agood representative of the underlying fiber bundles’ orientation, especially71.4. Anatomical and Functional Brain Connectivityat regions with complex fiber geometries, e.g. where fiber bundles meet orcross. For instance, when there are two crossing fiber bundles with similarradii at a given location, the diffusion tensor has a planar shape. In such acase, we cannot draw conclusions about the underlying microstructure sincedifferent bundle configurations may lead to the same planar shape.It is possible to characterize the diffusion process without imposing adiffusion model, such as the Gaussian model assumed in the case of diffusiontensors. It is apparent that this approach is theoretically less prone to biasesimposed by an explicit diffusion model [38]. The relationship between the 3DPDF of diffusion, also known as the spin propagator, and measured dMRIdata is given by the theory of q-space [39]. Given the spin propagator,the marginal probability of diffusion in any given direction, i.e. orientationdistribution function (ODF), provides a reduced representation of the PDFand is of interest for mapping white matter architecture. ODF can be foundusing a linear radial projection of the PDF [40] or by doing the projectionin a cone of constant solid angle (CSA) [41]. It has been shown that usingthe projection in a cone of CSA, denoted as CSA-ODF, results in a sharperrepresentation of the diffusion process leading to clearer recovery of certainfiber bundles [41].Using ODFs, which are more sensitive in capturing intricate intravoxelfiber geometries than diffusion tensors, may not directly translate into highersensitivity in AC estimation since the tractography method used also needsto make use of this extra information. For example, if only the dominantpeak of the ODF is used in a greedy manner in reconstructing white matterbundles, there may be no performance improvement at regions where a largefiber bundle crosses a smaller bundle. This is due to the partial volume av-eraging effect of diffusion tensors, i.e. the diffusion tensor would point alongthe more dominant bundle (the direction of the dominant ODF peak) at sucha location. On the other hand, additionally using the secondary ODF peakmay enable identifying the fiber crossing. However, this approach is alsoprone to errors since it necessitates estimating the number of bundles in avoxel to determine the number of ODF peaks to consider. A notable methodthat bypasses this potential source of error is global tractography [42], wherethe decisions are made on a global level as opposed to the greedy nature ofstreamline tractography. Global tractography reconstructs all fiber tractssimultaneously, which enables fiber trajectories to be jointly considered indetermining the most possible fiber configuration. In this approach, shortfiber segments are bridged together to form the set of fiber tracts that bestexplains the measured dMRI data. This global approach is particularly ben-eficial in resolving cases where local diffusion characteristics are ambiguous81.4. Anatomical and Functional Brain Connectivitydue to crossing fibers. Streamline tractography tends to terminate track-ing at such locations, whereas global tractography exploits the geometry ofthe surrounding fibers by considering the whole-brain fiber configuration inaggregate to compensate for the lack of reliable local information.Both streamline and global tractography provide us with point estimates,i.e. the variability of the reconstructed tracts based on underlying dMRIdata is not quantified. Several so-called probabilistic tractography methodshave been proposed to ameliorate this limitation. Probabilistic tractographytechniques are typically based on either modelling noise parameters by someprobability distribution [43, 44], or bootstrap-based methods that capturethe uncertainty in the data [45]. However, these methods are generally verycomputationally expensive, limiting their use when whole-brain connectivityis to be estimated. On a regular workstation, probabilistic tractographyof one brain at a typical resolution may take almost one day [46]. Theintensive computation involved thus inhibits the adoption of probabilistictractography in clinical tasks. Steps in a typical dMRI experiment are shownin Fig. 1.3.A major problem with AC estimation using tractography is the falsenegative connections due to premature termination of reconstructed tractsat fiber crossings. Given the recent evidence that 60% to 90% of the whitematter voxels contain crossing fibers [47], this limitation has serious impli-cations in terms of the quality of AC estimates. An example tractographyresult depicting this problem is given in Fig. 1.4, where the reconstructedtracts passing through corpus callosum are shown in red, and those passingthrough internal capsule are shown in green. Dotted red and green arrowsindicate the expected path of the tracts which could not be reconstructed.Such false negatives significantly hinder reliable AC estimation, especiallywhen estimating interhemispheric connections.Quantifying Anatomical ConnectivityPerforming tractography on dMRI data results in a set of reconstructed fibertracts. Unfortunately, this reconstruction is inherently non-quantitative interms of the connectivity strength. Although probabilistic tractography ap-pears to be more quantitative, it is important to realize that it providesestimates of our confidence on the fiber trajectories, and not connectionstrength. As an intuitive example, consider a case where brain regions Aand B are connected through a relatively small fiber bundle that has noimmediate neighboring bundles. Also consider regions C and D that areconnected through a major fiber bundle that crosses another fiber bundle91.4. Anatomical and Functional Brain ConnectivityDiffusion weighted imagesFiber tractsPreprocessingDiffusion modelingTractographyTensorsORODFsFigure 1.3: Sample workflow of a dMRI experiment. During a dMRI ex-periment, DW images sensitized to different directions are acquired, whereintensity of each voxel reflects the rate of diffusion at that location. Infor-mation from multiple directions are combined to reconstruct the diffusiontensors or ODFs. Fiber tracts are then delineated by following tensors orODFs with coherent orientations. MR scanner image: Big MRI by Liz West,CC BY 2.0, Anatomical and Functional Brain ConnectivityFigure 1.4: Depiction of the false negative problem in AC estimation usingstreamline tractography on diffusion tensors. Tracts extending from corpuscallosum (red) and internal capsule (green) are shown. Dotted arrows in-dicate the expected path of the tracts extending from corpus callosum andinternal capsule, which tractography failed to reconstruct. This figure isreproduced from [48] with the permission of the publisher.between C and D. In this hypothetical case, probabilistic tractography mightassign higher confidence to the connection between A and B, even thoughC and D are connected through a major fiber bundle that provides a highercapacity for interaction. We need quantitative measures of this capacity(which is loosely termed as AC strength here) for many problems such asdevising biomarkers of tissue structure and analyzing the relationship be-tween anatomical and functional connectivity in the brain.Quantifying AC in the brain is most commonly based on quantifyingone or more aspects of the fiber tracts reconstructed using streamline orglobal tractography, though computationally expensive approaches basedon probabilistic tractography techniques were also explored [49]. The choiceof which property to measure and of how to map it into an AC metric are keyaspects affecting the AC estimates. Arguably the most common AC metricis the number of reconstructed fiber tracts between pairs of brain regions,commonly referred to as fiber count. A variant of this approach involves anormalization of the fiber count by the total volume of the region pairs they111.4. Anatomical and Functional Brain Connectivityconnect to account for the variable size of the brain regions [50]. Besidesmetrics based on fiber count, use of the total length of reconstructed fibertracts has also been suggested as a measure of AC [51] aiming to correct forthe fact that longer tracts may have larger accumulated error, leading tolower fiber counts. Another metric, the average fractional anisotropy (FA)along streamlines connecting regions, has also been proposed as a proxy foranatomical connectivity strength [52]. FA is a local measure that describeshow anisotropic the underlying diffusion process is, and takes on valuesbetween 0 (perfectly isotropic diffusion) and 1 (an infinite cylinder).It is important to acknowledge that all of the aforementioned AC met-rics are confounded by several factors, limiting their interpretability. First,tractography can only delineate bundles of fibers in the brain, and not in-dividual fibers. The term fiber count can thus be misleading. Indeed, usingthe term streamline count has been recently proposed as an alternative [53].Nonetheless, we use the term fiber count for easier interpretation and toconform to the jargon used in existing literature, with the understandingthat it is the streamlines that are actually being counted. Moreover, wenote that the number of fibers is dependent on the number of seeds usedfor tracking the fibers, the tractography method used, and several featuresof the pathway such as curvature, length and width [53]. Additionally, wehighlight that FA not only depends on the reliability of local diffusion in-formation, but also on a large number of modulating factors such as axonalordering, axonal density, amount of myelination, and increase in extracel-lular or intracellular water [53]. However, such confounding factors did notimpede the adoption of a variety of AC metrics, driven by a practical needfor quantifying the degree of connection between brain regions. To this end,reconciling the presence of confounding factors with the practical need forquantitative connectivity estimation calls for a detailed analysis to deter-mine which AC metric has the highest potential of being of practical use inmultimodal brain image analysis efforts.1.4.2 Estimating Functional ConnectivityFC is conventionally defined as the temporal dependence of neuronal activ-ity patterns of remote brain regions [54], and it describes how brain regionpairs co-activate (positive synchrony, negative synchrony or no synchrony).Such co-activation patterns exist both during rest and task performance[22]. FC studies have been gaining traction in the clinical realm with appli-cations such as identifying group differences between populations, devisingbiomarkers for obtaining diagnostic and prognostic information, and guid-121.4. Anatomical and Functional Brain Connectivityance of invasive and noninvasive treatments [55].Several metrics have been proposed to quantify the dependence betweentime courses of the BOLD signal. Popular metrics include Pearson’s cor-relation [54], partial correlation [56], coherence [57] and partial coherence[57], with Pearson’s correlation being the most widely used metric by far.Pearson’s correlation measures marginal (direct and indirect) dependencies,whereas partial correlation quantifies conditional (only direct) dependencies.Coherence and partial coherence are the spectral counterparts of correlationmetrics. Though not commonly used, there exist more complex metricssuch as those based on higher-order statistics and temporal lag [58]. In sim-ulation based studies, it has been shown that correlation metrics performsignificantly better than more complex metrics [58]. On real data, test-retestanalyses confirmed that correlation is a more stable metric than coherence[59]. We have thus used Pearson’s correlation in this thesis to quantify FC.It is important to acknowledge that correlation between voxel time coursesis not informative of causality, or whether the connection is direct. For ex-ample, if brain region A is feeding into regions B and C, FC analysis basedon correlation would conclude that B and C are connected. Although thisbehaviour is consistent with the conceptual definition of FC given above(neural synchrony, or co-activation), it is important to understand the pre-cise nature of connectivity measured by the metric of choice. Steps in atypical fMRI experiment are shown in Fig. 1.5.A major problem with FC estimation is the many false positive connec-tions arising from confounds such as head movements, cardiac and respira-tory pulsations, and scanner drifts. Among these, head motion artifacts arelikely the most severe problem in fMRI studies [60]. Head motion resultsin artificial brain activations if it is not properly corrected for, since dataanalyses assume that each voxel represents a unique location in the brain.Such effects are typically visible as “ring activations” around the edges ofthe brain. Also, the voxels at brain boundaries may appear to be function-ally connected due to the factitious peak in brain activity they commonlyexhibit.1.4.3 Relationship between Anatomical and FunctionalBrain ConnectivityMotivated by how anatomical wiring of the brain facilitates and shapesfunctional interactions between brain regions, a growing interest in under-standing the relationship between AC and FC has emerged in recent years.To this end, AC-FC comparison has been done on a single axial slice [61],131.4. Anatomical and Functional Brain ConnectivityBOLDimagesFunctionalconnectivitymapVoxeltimecoursesPreprocessingConnectivity analysisPlot time seriesFigure 1.5: Sample workflow of an fMRI experiment. During an fMRIexperiment, the subject is either not engaged in any cognitive task (rest-ing state fMRI), or performs a series of cognitive tasks (task fMRI)while BOLD images are being acquired. The preprocessed voxel timecourses are compared with each other to infer the functional connectiv-ity structure. MR scanner image: Big MRI by Liz West, CC BY 2.0, Functional Segregation in the Brain: Activation Detectionwithin a specific subnetwork [8, 62], and for the whole brain [50, 63, 64]. Col-lectively, these studies suggest that there is a robust and strong correlationbetween estimates of AC and FC. Although these earlier works comparingAC and FC suggest high consistency between their estimates, it has beenmore recently shown that the AC-FC correlation is rather low with typicaldata acquisitions [65]. A fundamental factor governing the level of consis-tency is the way in which FC and AC are estimated. Although neuroimagingcommunity largely settled on the use of Pearson’s correlation between voxeltime courses for estimating FC [58, 59], investigations of AC estimation andquantification have received less attention. False negatives are common intractography mainly due to the abundance of crossing fibers in the humanbrain. There are two potential solutions to ameliorate this problem. Thefirst one involves increasing the spatial resolution of dMRI data such thateach voxel covers a smaller region in the brain, reducing the probability thata given voxel will comprise multiple fiber bundles. The second solution isto employ a tractography algorithm that handles fiber crossings better thanthe widely used streamline tractography. Furthermore, it is important todetermine which AC metric has the highest potential of being of practicaluse in multimodal brain image analysis efforts.1.5 Functional Segregation in the Brain:Activation Detection1.5.1 Unimodal ApproachesThe aim of activation detection is to map brain regions to function by inves-tigating the regions that are activated in response to stimuli. The standardapproach for inferring brain activation from fMRI data is to model the fMRIobservations at each voxel independently as a linear combination of expectedtemporal responses using the general linear model (GLM) [66] to estimatethe voxel-wise activation effects. Statistics, such as t-values, correspondingto the p-values of observing the estimated activation effects when there is noactivation are combined in an activation statistics map. This map is thenthresholded, resulting in a set of supra-threshold voxels that are declaredas active. However, since inference is done for all voxels separately, thereis a need for adjusting the threshold for multiple comparisons [67]. A ma-jor shortcoming of the GLM approach is that, spatial correlations amongvoxels are ignored by independently testing each voxel. To account for thetendency of the nearby voxels to work in tandem, activation statistics maps151.5. Functional Segregation in the Brain: Activation Detectionare typically smoothed with a Gaussian kernel as a post-processing step[68]. Bayesian schemes incorporating spatial priors to smooth the data ina data-driven way have also been proposed [69, 70]. These adaptive ap-proaches have been shown to preserve the spatial features of the activationpatterns that would be lost with a fixed Gaussian kernel [70]. Evidently,both adaptive and nonadaptive filtering approaches result in reduced spatialresolution, decreasing the localization power of the analyses. Alternatively,Markov random fields (MRF) based approaches which view activation de-tection as a class labelling problem, i.e. active and nonactive, have also beenexplored [71–73]. These methods encourage spatial continuity by penalizingthe discrepancies of labels among neighbouring voxels, as opposed to refin-ing activation statistics maps. Although such methods help suppress falsespatially-isolated activations, they completely ignore long-range functionalinteractions. To address this limitation, incorporation of functional connec-tivity information into activation detection has been put forth [74], but theapproach taken estimates both activation effects and functional connectiv-ity from the same dataset, hence the information gain might be limited.Given the typically strong noise in fMRI data, exploring other sources ofinformation to regularize activation detection may be beneficial.1.5.2 Multimodal ApproachesA major challenge in brain activation detection is the existence of sponta-neous fluctuations in the brain in addition to task-evoked signal [75]. Suchspontaneous modulations in brain activity are continually present at alltimes [76]. Thus, part of the noise observed in task fMRI data ascribesto ongoing RS activity. As opposed to the traditional view that RS activityis just noise that needs to be averaged out, strong synchrony in RS modu-lations between specific brain areas has been observed in numerous studies[77]. In fact, many of the detected RS subnetworks exhibit high resem-blance to subnetworks seen in task experiments [78]. Further supportingthis finding is a recent work [79] that demonstrated enhanced sensitivity intask activation detection by incorporating an RS connectivity prior. Takentogether, these findings suggest that there is spatial structure in RS activitywhich is likely similar at rest and during task performance [75].Since the time at which RS activity peaks and troughs is difficult topredict, decoupling task-evoked activity from RS fluctuations is non-trivial.To the best of our knowledge, the only previous work that attempted totackle this challenging problem of RS activity removal from task fMRI datawas by Fox et al. [75]. Specifically, the authors showed that for a right161.6. Functional Integration in the Brain: Subnetwork Identificationhanded motor task, subtracting out fMRI signals in the right somatomotorcortex (RSC) from the left somatomotor cortex (LSC) after weighting by aregression coefficient derived from the RS data significantly reduced inter-trial variability in fMRI response, suggesting that RS fluctuations and task-evoked signal are superimposed roughly in a linear fashion in the humanbrain. The basis of this approach is twofold. First, the presence of coherentRS activity between LSC and RSC is well established [22]. Second, righthanded motor tasks typically activate only LSC, thus signals in RSC wouldlargely correspond to RS activity. Subtraction of signals in RSC from LSCwould thereby remove the RS components within the task fMRI time coursesof LSC. However, not all tasks evoke only lateralized activation. Thus,simply subtracting signals in one side of the brain from the other side is notalways suitable for removing RS modulations from task fMRI data.1.6 Functional Integration in the Brain:Subnetwork Identification1.6.1 Unimodal ApproachesSince the interactions between specialized regions in the brain are knownto play a key role in mediating brain function [80], fMRI has been widelyused for the modelling of voxel interactions. One of the main approachesfor studying this integrative property of the brain is to group voxels intofunctional subnetworks based on the temporal dependencies among theirfMRI observations. The most widely used techniques for this include theseed-based approach and independent component analysis (ICA) [81]. Inthe seed-based approach, functional subnetworks are estimated based onthe temporal correlations between a seed region and all other regions in thebrain [81]. Since pre-specifying a seed might be difficult in some cases, tech-niques such as k-means [82], hierarchical clustering [83], and Gaussian mix-ture model (GMM) [84] have been proposed to automatically identify seedregions. Despite its simplicity and intuitiveness, the seed-based approachhas the drawback that results are restricted by the choice of seeds [85]. Incontrast, ICA provides a data-driven means for extracting subnetworks bydecomposing fMRI observations into maximally independent spatial com-ponents and associated time courses [86]. At the intra-subject level, voxelsassigned weights larger than a threshold within the same spatial componentare grouped as a subnetwork. However, the threshold is generally chosen inan ad hoc manner, which precludes statistical interpretation [86]. In prac-171.6. Functional Integration in the Brain: Subnetwork Identificationtice, ICA is typically performed at the group level by pooling informationfrom multiple subjects, in which a statistically meaningful threshold can beset [87].An important point commonly overlooked in the literature is the factthat certain brain regions are known to interact with multiple subnetworks[88]. It is thus important to develop methods that allow spatial overlapamong identified subnetworks. Recently, a number of techniques have beenput forth for identifying overlapping brain subnetworks. In [89], the tech-nique involves first performing maximal clique detection and treating eachmaximal clique as a node of a new graph. Subnetworks are then extractedby applying modularity optimization [90] on this new graph. Since the max-imal cliques might share nodes, overlaps between subnetworks are enabled.The main shortcoming of this technique is that it requires graph edges to bebinarized since it cannot handle weighted graphs, but brain connectivity lieson a continuum. This need for edge binarization could thus result in impor-tant information being discarded. A technique adopted from social networkanalysis called connected iterative scan (CIS) has also been explored forbrain subnetwork identification [91]. Taking each node as a candidate sub-network, CIS determines if other nodes belong to this candidate subnetworkbased on their influence on a graph density metric. A limitation of CIS isthat it is sensitive to a weighting factor that controls subnetwork size. An-other technique for finding overlapping brain subnetworks is temporal ICA(tICA), which is more suited for finding overlaps than the more widely-usedspatial ICA (sICA), since maximizing spatial independence tends to resultin little overlaps [88]. However, the proposed implementation uses an ad hocthreshold for extracting interacting nodes from the spatial maps, which lim-its the statistical interpretability of the results. Grouping the edges betweennodes instead of the nodes themselves has also been proposed for detectingsubnetwork overlaps [92]. In that technique, edges are grouped using e.g.hierarchical clustering with the set of nodes associated with any edge of agiven edge cluster declared as a subnetwork [93]. The main limitation ofthis technique is scalability since in weighted graphs where all nodes areconnected by an edge, the number of pairwise relationships between edgesgrows quartically with the number of nodes, e.g. 100 nodes would have morethan 108 relationships between edges. Several other overlapping communitydetection techniques have been proposed in the general context, a review ofwhich can be found in [94]. However, to the best of our knowledge, noneof these overlapping subnetwork identification methods permit multimodalintegration.Another limitation in most community detection techniques is their sole181.6. Functional Integration in the Brain: Subnetwork Identificationfocus on optimizing a criterion, such as modularity, without considering thevariability of the subnetwork assignment under data perturbation or whatsubnetworks one would get from random graphs to assess statistical signifi-cance. In the area of brain activation detection, rigorous statistics has longbeen employed to control for false positives [67]. It is known that even ran-dom networks can have high modularity [95], thus statistically controllingfor false node inclusion in subnetwork identification is crucial. A notableexception is the order statistics local optimization method (OSLOM) [96].OSLOM uses a statistical criterion based on a null model derived from ran-dom graphs for assessing the significance of the identified subnetworks, and isshown to outperform many state-of-the-art community detection techniques[94].1.6.2 Multimodal ApproachesIncorporating anatomical information extracted from dMRI data into theinvestigation of functional brain dynamics has attracted growing interestsince functional synchronization between spatially distinct brain regions isenabled through neural fiber pathways [97, 98]. Most of the early workfocused on direct comparisons of anatomical and functional connectivitylearned separately from dMRI and fMRI data as explained in Section 1.4.3.The main findings were that brain areas with high AC typically displayhigh FC but the converse does not always hold [98]. However, recent studieshave shown that the correlation between FC and AC estimates is ratherlow with typical data acquisition schemes [65]. As explained in Section1.4, FC estimates tend to contain many false positive connections whereasAC estimates often suffer from false negatives. These discrepancies mayappear as setbacks at first sight, but these inter-modality differences in factrender multimodal integration useful since incorporating AC would reducethe effects of noise-induced FC and vice versa.The inherent relationship between anatomy and function in the brainrevealed by the aforementioned studies implicates the potential benefits offMRI-dMRI fusion such as improving the sensitivity of data analyses usingcomplementary information or improving our understanding of how brainfunction is facilitated by anatomy in general. In line with this observation,multimodal integration for joint anatomical and functional connectivity in-ference have been recently explored to better delineate the architecture ofthe human brain [99, 100]. It has been shown that multimodal connectivityinference offers higher inter-subject consistency of the detected connectionpatterns [99] and higher discrimination power between the control and clin-191.7. Thesis Contributionsical populations [100] compared to analyzing each modality in isolation.In addition to multimodal connectivity inference studies, multimodalsubnetwork extraction problem has also recently started to attract atten-tion. For this problem, a multi-view spectral clustering algorithm has beenrecently adopted in [52], where FC and AC are considered as two views ofthe brain. The Laplacian of the similarity matrix of one view is projectedonto the eigenspace of the other view in an alternating manner, with spec-tral clustering subsequently applied on the converged eigenspace to identifysubnetworks. A related approach has been used in [101], where the commoneigenspace between the similarity matrices of FC and AC is found via jointdiagonalization, on which spectral clustering is performed.1.7 Thesis ContributionsWe propose in this thesis novel directions for multimodal brain image anal-ysis for assessing functional segregation and integration in the human brain.We demonstrate the strengths of multimodal approaches over traditionalunimodal analyses on carefully designed synthetic datasets as well as realfMRI and dMRI datasets. Our contributions in advancing the state-of-the-art are listed below. We note that the original mathematical notationsused in publications that led to this thesis are kept to ease cross-referencing.Hence, there might be some overlaps in notations between different sections.To resolve any ambiguities that might arise due to such overlaps, we haveexplicitly redefined the notations in each subsection.1.7.1 Anatomical and Functional Brain Connectivity• We showed that the consistency between estimates of AC and FC canbe increased using streamline tractography on ODFs as opposed toconventional diffusion tensors when estimating AC. We also demon-strated a further improvement when global tractography on ODFs isdeployed [P1].• We showed that part of the inconsistency between estimates of AC andFC can be attributed to the relatively low spatial resolution of dMRIdata, often leading to partial volume effects. We proposed to use asuper-resolution approach based on dictionary learning for alleviatingthis problem and showed that our method outperforms traditionalinterpolation strategies [P2].201.7. Thesis Contributions• We showed that the way AC is quantified has a significant effect on theagreement between estimates of AC and FC. We demonstrated thatfiber count and volume-normalized fiber count are the two metrics thatcorrelate best with FC, and that fiber count has lower bias in favor ofshorter fibers than volume-normalized fiber count [P3]. Fiber count isthus used to quantify AC in this thesis.1.7.2 Multimodal Brain Activation Detection• We proposed a novel approach that exploits the intrinsic task-restrelationships in brain activity to reduce the confounding effects of RSactivity on task activation detection. We showed that removal of theestimated RS modulations from task fMRI data significantly improvesactivation detection [P4].• We proposed incorporating AC into brain activation detection to mit-igate the effects of noise. We illustrated that incorporating AC sig-nificantly increases sensitivity in detecting task activation for healthysubjects. We also showed that this multimodal approach enables thedetection of significant group differences between schizophrenia pa-tients and controls that are missed with standard methods [P5].1.7.3 Multimodal Brain Subnetwork Identification• We proposed a multimodal integration technique for identifying sub-networks of brain regions that exhibit high inter-connectivity bothfunctionally and anatomically. Our proposed technique allows forsubnetwork overlaps, and provides statistical control on false nodeinclusion in the identified subnetworks. We demonstrated that ourtechnique achieves significantly higher subnetwork identification accu-racy than state-of-the-art techniques on synthetic data. On real data,we showed that our method attains improved test-retest reliability onmultiple network measures. Further, we showed that the tasks subjectsare involved in can be more accurately classified based on the subnet-works employed during different tasks with this multimodal approachcompared to unimodal approaches [P6-P8].21Chapter 2Anatomical and FunctionalBrain ConnectivitySince anatomical wiring of the brain shapes functional interactions, it is intu-itive to expect a strong correlation between AC and FC in the brain. Indeed,a number of studies suggested a strong positive correlation between AC andFC estimates [63, 64, 98]. However, it has been more recently shown thatthe AC-FC correlation is as low as ∼0.12 using the most common metricsto quantify AC and FC (fiber count for AC, Pearson’s correlation betweenfMRI time courses for FC) on typical dMRI (∼32 gradient directions) andRS fMRI (∼7 min) data [65]. We hypothesize that this result is in partdue to error propagation during the several steps of AC and FC estimation.As explained in Section 1.4.2, the neuroimaging community has largely set-tled on the estimation of FC using Pearson’s correlation between fMRI timecourses following a series of standard preprocessing steps. However, AC es-timation using dMRI data is far from being standardized. We thus focusedour attention on various steps leading to quantitative AC estimates.Firstly, accuracy of the reconstructed fiber tracts is often hampered bythe inherently low resolution of dMRI data. Currently achievable spatialdMRI resolution is around 2×2×2 mm3, while the actual neuronal fiber di-ameter is on the order of 1 µm [102]. A voxel can thus comprise severaldistinct fiber bundles with differing orientations, leading to partial volumeaveraging [103]. At such locations, diffusion information typically becomesambiguous, and tractography is often falsely terminated. Therefore, in-creasing the spatial resolution in dMRI holds great promise towards moreaccurate delineation of fiber tracts. To this end, we propose to use a super-resolution approach based on dictionary learning for alleviating this problem(Section 2.1). Unlike the majority of existing super-resolution algorithms,our proposed solution does not entail acquiring multiple scans from the samesubject which renders it practical in clinical settings and applicable to legacydata. Moreover, this approach can be used in conjunction with any diffusionmodel.Secondly, given a set of DW images, the diffusion model and tractogra-22Chapter 2. Anatomical and Functional Brain Connectivityphy algorithm employed directly affect the quality of the reconstructed fibertracts. Suboptimal choices made at this step may result in erroneous ACestimates regardless of the quality of the underlying data. dMRI data haveoverwhelmingly been modelled using tensors to date, which approximatesthe local diffusion process by a single Gaussian probability density function.Due to this simplistic approximation, diffusion tensors can only capture theprincipal diffusion direction, which makes them prone to errors induced bycrossing fibers. With recent advances in dMRI acquisition, a larger numberof gradient directions can now be acquired within a practical amount of time,which enables more sophisticated methods that can model multiple fiber di-rections, such as ODF reconstruction [15], to be readily employed. To betterexploit the high angular resolution data, ample efforts have been placed onimproving ODF estimation strategies during recent years. However, less at-tention has been paid to tractography methods, which are equally importantfor AC estimation. Streamline tractography performed on diffusion tensorsor ODFs remains to be the norm. Typically, streamline tractography fol-lows either the principal diffusion directions of diffusion tensors or largestODF peaks in tracking fibers one at a time. Approaches that consider allfibers in aggregate have recently been explored [42]. Whether these globaltractography techniques can help increase the consistency between AC andFC estimates has to date not been investigated. To this end, we analyze theimpact of AC estimation strategies on the consistency between AC and FCestimates. Specifically, we show increased AC-FC consistency with stream-line tractography on ODFs compared to using conventional diffusion tensors.We also demonstrate a further improvement when global tractography onODFs is deployed (Section 2.2).Lastly, the choice of which streamline property to measure and of how tomap it into an anatomical connectivity metric are key aspects affecting theAC estimates once the streamlines are reconstructed. To determine whichAC metric has the highest potential of being of practical use in multimodalbrain image analysis efforts, we compare four commonly used AC metrics interms of their impact on the relationship between estimates of AC and FC.We demonstrate that fiber count and volume-normalized fiber count are thetwo metrics that correlate best with FC, and that fiber count has lower biasin favor of shorter distances than volume-normalized fiber count (Section2.3). Fiber count is thus used to quantify AC in this thesis.232.1. Super-resolution for dMRI2.1 Super-resolution for dMRIDue to partial volume averaging, local dMRI observations typically becomeambiguous in voxels comprising several distinct fiber bundles, and tractogra-phy is often falsely terminated. Therefore, increasing the spatial resolutionin dMRI holds great promise towards more accurate delineation of fibertracts. There are practical limitations in increasing the resolution of theacquired data directly, such as reduced SNR and prolonged scanning time[104]. Such limitations motivate the search for post-processing solutions forincreasing the spatial resolution, such as super-resolution techniques.Super-resolution techniques have been previously adopted to increasethe spatial detail in dMRI. In the literature, the term super-resolution isused for two distinct classes of methods which follow different paradigms.The first class of methods are based on performing multiple low-resolutionacquisitions, followed by the fusion of the information in these images togenerate high-resolution images. To this end, fusing images spatially shiftedat subvoxel level [105], as well as fusing multiple anisotropic images withhigh resolution only along one axis [102, 106] have been explored. In a fairlysimilar spirit, combining DW images acquired at two different resolutions toinfer high-resolution diffusion parameters using a Bayesian model has alsobeen proposed [107]. The inherent drawback of these approaches is the de-pendence on a specific acquisition protocol, limiting their usability in generalsettings. The second class of methods do not require multiple acquisitions,and these are typically based on examples or priors about the correspon-dence between low and high resolution images. Falling in this category, anapproach to reconstruct diffusion tensors at a resolution higher than theunderlying DW images using a single dMRI acquisition has been recentlyproposed [108]. Even though this method eliminates the need for multipleacquisitions, it is only geared towards estimating diffusion tensors, and can-not be easily extended to higher order diffusion models such as ODFs. To thebest of our knowledge, the only previous work that tackled the problem ofsuper-resolving dMRI data from a single acquisition independent of the dif-fusion model was by Coupe´ et al. [109]. Specifically, the authors showed thatsuper-resolving b=0 (non-diffusion-weighted) image using a locally adaptivepatch-based strategy, and using this high-resolution b=0 image to drive thereconstruction of DW images outperforms upsampling of dMRI data usingclassical interpolation methods. Beyond these two classes, a new perspec-Section 2.1 is adapted from: B. Yoldemir, M. Bajammal, R. Abugharbieh. Dictionarybased super-resolution for diffusion MRI. In Proceedings of the MICCAI Workshop onComputational Diffusion MRI, pp. 194-204, Boston, MA, September 2014.242.1. Super-resolution for dMRItive to gain spatial resolution in dMRI has been proposed which is termed assuper-resolution tract-density imaging (TDI) [110]. This approach is funda-mentally different than the aforementioned super-resolution methods in thesense that the aim is to generate high resolution tract density maps throughcounting the number of tracts present in each element of a subvoxel grid,rather than super-resolving the DW volumes prior to tractography.We start by presenting our assumed data acquisition model (Section2.1.1). Given a set of acquired DW volumes, we form a training set that in-cludes the original volumes and a set of corresponding downsampled volumesat double the voxel size. We then construct two over-complete dictionariesfrom the original and downsampled set of volumes (Section 2.1.2). For apreviously unseen input DW volume, we obtain the super-resolution data intwo steps. First, we sparsely code the volume against the dictionary learnedfrom the downsampled volumes in the training set. We finally apply thegenerated sparse code to the dictionary learned from the original resolutionset to obtain the super-resolution data (Section 2.1.3).2.1.1 Acquisition ModelLet vL be an acquired volume and vH be the corresponding unobservedhigher resolution volume. We assume that the relationship between thesetwo volumes is modelled by [111]vL = S B vH + n, (2.1)where S is a downsampling operator, B is a blurring operator and n isadditive white Gaussian noise. We aim to invert this acquisition modelto approximate the unobserved higher resolution volume through super-resolution. The maximum-likelihood solution to this problem involves theminimization of ||S B vˆH − vL||2, where vˆH is the estimated high resolutionvolume. However, the inversion of S B is ill-posed [111], hence infinitelymany maximum-likelihood solutions exist. We thus cast the problem in adictionary learning framework instead, as explained next.2.1.2 Dictionary ConstructionWe model each 3D patch in dMRI volumes as a sparse linear combination ofatoms from a learned dictionary D. In the proposed approach, we use twodictionaries to capture the correspondence between low and high resolutiondMRI volumes. These two dictionaries are learned from the original trainingdataset and its downsampled version, respectively.252.1. Super-resolution for dMRILet vO be the set of original training volumes concatenated across scansand vD be the corresponding set of downsampled volumes. We extract alloverlapping patches in these two sets of volumes, denoted by pO and pD,respectively. Using pO and pD, we construct two over-complete dictionariesas follows:minDO,DD,y∑‖pD −DDy‖22 +∑‖pO −DOy‖22 + ψ(y), (2.2)where y = {y(i,j,k)} is the set of sparse coding vectors for each image loca-tion (i, j, k), and DD and DO are the generated over-complete dictionariesof the downsampled and original volumes, respectively [111]. ψ(y) is a reg-ularization term which we set to be ψ(y) = ||y||1, inducing sparsity on thegenerated coding vector [112]. We note that the same set of coding vectorsy is used for both dictionaries. In other words, the learned atoms of thetwo dictionaries represent matched pairs. We set the number of atoms ineach dictionary to 1000 and the patch size to 3×3×3 voxels, which wereempirically chosen to strike a balance between representation accuracy andoverfitting. Note that these two parameters can be optimized via cross-validation, however we did not perform such analysis since our aim is tohighlight that AC-FC agreement is affected by the relatively low spatial res-olution of dMRI data rather than achieving the highest correlation betweenthe estimates of these two types of connectivity.2.1.3 Super-Resolved Volume GenerationLet pI be the set of low resolution overlapping patches obtained from apreviously unseen input volume that we wish to super-resolve. We code pIwith respect to DD asminyI||pI −DDyI ||22 + ψ(yI), (2.3)where yI is the set of coding vectors for pI , with ψ(yI) again being the l1norm of yI , enforcing sparsity on the coefficients. Once the input volume issparsely coded using DD, we generate a new set of super-resolved patches,pS , by applying the sparse coding vector yI to DO previously constructedfrom the training data:pS = DOyI . (2.4)We note that this process results in a patch being generated for eachvoxel. We then reconstruct the super-resolved volume by averaging neigh-bouring overlapping patches. We used K-singular value decomposition (K-SVD) [113] to construct the dictionaries and orthogonal matching pursuit262.1. Super-resolution for dMRI[114] to sparsely code the 3D patches. Theoretically, pO, pD and pI canbe extracted at once from the volumes of all gradient directions in the DWimages. However, we opt to apply super-resolution for each gradient direc-tion separately. This helps circumvent computational limitations that mightarise, especially with the increasingly large number of gradient directions ac-quired in practice.2.1.4 ExperimentsMaterialsWe validated our method on the publicly available multimodal Kirby 21dataset.1 Along with other imaging modalities, this dataset comprises dMRIand RS fMRI scans of 21 subjects with no history of neurological disease(11 men, 10 women, 32±9.4 years old). We summarize the key acquisitionparameters below. Further details on data acquisition can be found in [115].In our experiments, we used 10 subjects for dictionary training, and 10 othersubjects for testing.RS fMRI data of 7-minute duration were collected with a TR of 2 s and avoxel size of 3 mm (isotropic). Preprocessing steps included motion correc-tion, bandpass filtering at 0.01 and 0.1 Hz, and removal of white matter andcerebrospinal fluid confounds. A common approach to brain image analysisinvolves examining regions of interest (ROIs) as internally coherent units.The representative signal of a given ROI is defined either as the average timecourse of the voxels in the ROI or the first principal component from a prin-cipal component analysis (PCA) performed on the time courses of individualvoxels [116]. The main reasons for working on voxel groups instead of vox-els include reducing the data dimensionality for computational tractabilityand controlling for type-I error by reducing the number of statistical tests[117]. Since each voxel is unlikely to function in isolation [118], ROI-basedapproach offers an efficient alternative to voxel-based analyses. We thusadopted an ROI-based approach instead of a voxel-based approach in thisthesis. We divided the brain into 150 parcels by applying Ward clustering[119] on the voxel time courses, which were temporally concatenated acrosssubjects. Parcel time courses were then found by averaging the voxel timecourses within each parcel.The dMRI data had 32 diffusion-weighted images with a b-value of 700s/mm2 in addition to a single b=0 image, with a voxel size of 0.83×0.83×2.2mm3. Since anisotropic voxels were previously shown to be suboptimal for1This dataset is available online at: Super-resolution for dMRIfiber tractography [120], we resampled each volume to 2 mm isotropic res-olution prior to any analysis. We also applied a Rician-adapted denoisingfilter [121] to eliminate nonstationary noise commonly observed in DW im-ages, since our acquisition model assumes Gaussian noise. We then warpedour functionally derived group parcellation map to the b=0 volume of eachsubject using FSL [31] to facilitate the computation of fiber count.ResultsWe first present a qualitative comparison between the fiber tracts recon-structed from the original (2 mm isotropic) and super-resolved (1 mm isotrop-ic) dMRI data. For ease of interpretation, we chose to employ deterministicstreamline tractography with the diffusion tensor model, which is by farthe most popular tractography approach to date. However, we highlightthat our super-resolution approach can be used with any diffusion modeland any tractography method. Tractography was carried out using Dipy[122], with 750,000 seed points for both the original and super-resolutiondata. We generated the tract-density maps by calculating the total num-ber of fiber tracts present in each voxel. Fig. 2.1(a),(c) and (b),(d) showsample tract-density maps with the original and super-resolved dMRI data,respectively. As observed from these figures, the tract-density maps gener-ated from the super-resolution data clearly convey more spatial information.Fig. 2.1(e),(f) and (g),(h) show the corticospinal tracts extracted using anROI placed on the brain stem for two representative subjects. It can beobserved that fiber tracts reconstructed from the super-resolution data cancapture the fan-shape configuration of the corticospinal tract more fully.To quantify the improvement in tractography with the suggested ap-proach, we analyzed the consistency between measures of intra-subject ACand FC. We estimated AC using the fiber counts between brain region pairs,and FC using Pearson’s correlation between parcel time courses. For eachsubject, AC and FC are vectors of size d(d − 1)/2 × 1 comprising the cor-responding connectivity estimates between each region pair, where d is thenumber of brain regions. We then calculated Pearson’s correlation betweenintra-subject AC and FC to quantify the consistency between the two con-nectivity estimates. This linear model was used instead of higher-ordernonlinear models in estimating AC-FC consistency for easier interpretabil-ity, and based on the empirical finding that AC and FC in the human brainare roughly linearly related [64, 97]. Using this correlation measure, wecompared the proposed super-resolution approach with trilinear and splineinterpolation in addition to an alternative super-resolution method; collab-282.1. Super-resolution for dMRIOriginal	   Super-­‐resolved	  (a)	   (b)	  (c)	   (d)	  (e)	   (f)	  (g)	   (h)	  Figure 2.1: Qualitative comparison between the tract-density maps andfiber tracts reconstructed from the original (left) and super-resolved (right)dMRI data. Original dataset has 2 mm isotropic resolution which is super-resolved to 1 mm isotropic resolution. Each row corresponds to a differenttest subject. Tract-density maps of super-resolved data ((b) and (d)) showmarkedly improved spatial detail compared to those of original data ((a)and (c)). Corticospinal tracts reconstructed from super-resolved data ((f)and (h)) can capture the fan-shape configuration more accurately than thosegenerated from original data ((e) and (g)).292.1. Super-resolution for dMRIorative and locally adaptive super-resolution (CLASR) [109]. To the bestof our knowledge, CLASR is the only existing single image super-resolutionmethod for dMRI which is independent of the diffusion model employed.Fig. 2.2 shows the AC-FC correlation for each subject tested. Taking theaverage AC-FC correlation across the group when using the original dataas a baseline, the improvement was 5.7% with spline interpolation, 13.6%with CLASR, and 27.1% with our proposed method. On the other hand,there was a 6.3% decrease in the correlation when trilinear interpolationwas used. The difference in the performance of our method and every othermethod tested was found to be statistically significant at p<0.01 based onthe Wilcoxon signed-rank test, showing its potential for enhanced anatomi-cal connectivity assessment. Our results thus suggest that low spatial resolu-tion of dMRI data can partially account for the low AC-FC correlation, andstatistically significant improvements can be achieved using super-resolveddMRI data.To investigate why trilinear interpolation resulted in a lower AC-FCcorrelation compared to the original data, we calculated the number of tractsreconstructed with each method. The local intra-parcel connections wereexcluded since they have no effect on AC-FC correlation. Fig. 2.3 showsthe number of inter-parcel tracts averaged across the group along with thecorresponding standard deviations. As observed from this figure, performingtractography on volumes upsampled with trilinear interpolation resulted ina lower number of tracts compared to the original volumes, even thoughthe same number of seed points were used to initiate tracking for all of themethods we compared. We speculate that the reason of this phenomenonis the additional partial volume effects introduced by the blurring of thedata during trilinear interpolation, which hamper the tractography quality.Spline interpolation, however, is known to cause less blurring compared totrilinear interpolation, and our results suggest that upsampling dMRI datausing spline interpolation can be beneficial for tractography. The overalltrend of inter-parcel tract counts closely resembles to that of the AC-FCcorrelation, with our proposed method outperforming all other methodstested. This shows that dictionary based super-resolution is a viable post-processing solution for dMRI that can help in mapping the white matterbrain architecture more accurately.2.1.5 DiscussionLow spatial resolution is a known limitation of dMRI, which often hindersthe performance of tractography significantly. We proposed the use of a302.1. Super-resolution for dMRI1 2 3 4 5 6 7 8 9  Student Version of MATLABOriginal Trilinear Spline CLASR Our methodAC−FC CorrelationSubject IDFigure 2.2: AC-FC correlation for 10 subjects with AC estimated from thedata at its original resolution (2 mm isotropic), and high-resolution data(1 mm isotropic) obtained using trilinear interpolation, spline interpolation,CLASR and the proposed method. Our method outperforms all other meth-ods tested for eight of the subjects, and performs comparable to CLASR fortwo subjects (subjects 4 and 10).312.1. Super-resolution for dMRI4.555.566.57x 104Student Version of MATLABInter−parcel tract countOriginal Trilinear Spline   CLASR Our methodFigure 2.3: Number of inter-parcel tracts reconstr cted from the dataat its original resolution (2 mm isotropic), and high-resolution data (1mm isotropic) obtained using trilinear interpolation, spline interpolation,CLASR and the proposed method. Intra-parcel tracts are not included heresince they do not contribute to AC-FC correlation. We emphasize thattractography is initiated with the same number of seeds (750,000) for eachmethod.322.2. Estimating Anatomical Connectivitysimple yet effective super-resolution technique in dMRI to capture a moreaccurate portrayal of the white matter architecture. The advantage of thismethod is three-fold. First, the super-resolved DW images can be used withany diffusion model as permitted by the number of gradient directions inthe original dataset. Second, this method does not rely on repeated acqui-sitions from the same subject, allowing it to be used with legacy data andunder various clinical acquisition schemes. Third, this method may still bereadily applied when the imaging protocol involves multiple acquisitions, asan additional step after reconstructing a single image from multiple low res-olution acquisitions [102, 105, 106]. Quantitatively, we demonstrated thatAC-FC consistency can be markedly increased with the use of our approachin estimating AC. We also qualitatively illustrated that the gain in spatialresolution remarkably improves the fiber tracts and tract-density maps gen-erated. Taken collectively, our results suggest that dictionary based super-resolution holds great promise in enhancing the spatial resolution in dMRI,without requiring additional scans or any modifications of the acquisitionprotocol.2.2 Estimating Anatomical ConnectivityThe success of techniques that integrate dMRI and fMRI for multimodalconnectivity estimation highly depends on the consistency between the ACand FC estimates. To this end, we analyze the impact of AC estimationstrategies on the consistency between AC and FC estimates. Specifically,we compare conventional streamline tractography on conventional diffusiontensors and ODFs. Our rationale is that ODFs can model multiple fiber di-rections, thus may lead to increased AC-FC consistency by more accuratelydepicting the anatomical wiring structure of the brain. Further, we comparestreamline tractography against global tractography [42]. Global tractogra-phy accounts for spatial interactions among fiber orientation estimates byreconstructing all fibers simultaneously and finding the configuration thatbest explains the data [42]. Global tractography thus alleviates error ac-cumulation along reconstructed tracts, which is a major drawback of thestreamline approach. We hence hypothesize that fiber tracts extracted withglobal tractography would better resemble the macroscale neural circuitry,Section 2.2 is adapted from: B. Yoldemir, B. Ng, R. Abugharbieh. Effects of tractog-raphy approach on consistency between anatomical and functional connectivity estimates.In Proceedings of the IEEE International Symposium on Biomedical Imaging (ISBI’14),pp. 250-253, Beijing, China, April 2014.332.2. Estimating Anatomical Connectivityand thus provide higher consistency between AC and FC estimates. We usethe most common methods for estimating FC and quantifying AC, namelyPearson’s correlation between fMRI time courses and fiber count, respec-tively. We focus our attention on AC estimation, i.e. diffusion modellingand tractography.2.2.1 Modelling of Local Diffusion PropertiesWhen local diffusion is modelled using diffusion tensors, only the dominantdiffusion direction can be recovered at each voxel. In contrast, ODFs aremore sensitive in capturing intricate intravoxel fiber geometries. ODFs aretypically estimated from Q-ball imaging reconstruction of high angular reso-lution diffusion imaging (HARDI) data acquired on one or multiple sphericalshells in q-space [15]. We opt to use CSA-ODF [41] since this reconstructionmethod intrinsically provides sharp ODFs, which enables multiple intravoxelfiber orientations to be more easily resolved. Moreover, this method takesadvantage of the additional information provided by multi-shell acquisition,which has been empirically shown to better resolve fiber crossings comparedto single q-shell models [41]. Using CSA-ODF thus better exploits the mul-tiple q-shells available in the data used in our experiments. Further detailson CSA-ODF estimation can be found in [41].2.2.2 TractographyStreamline TractographyConventional streamline tractography is a deterministic approach, in whicha single fiber trajectory is generated for each seed point [36]. Seed pointsare defined on a dense grid at the voxel or subvoxel level. Fiber tracts arethen reconstructed by bidirectionally traversing the direction of maximumdiffusion, i.e. the largest ODF peak or the principal diffusion directionin the case of diffusion tensors. In this work, we use Euler’s method fornumerically performing tractography. Starting from a seed point p0, fibertracts are reconstructed usingpn+1 = pn + v(pn)∆s, (2.5)where pn is the position of the propagator at step n, v(pn) is the propagationdirection at pn, and ∆s is the fixed step size [122]. To generate smoothtracts, we compute v(pn) using trilinear interpolation which combines theinformation of voxels within the 8-neighbourhood of pn [122].342.2. Estimating Anatomical ConnectivityGlobal TractographyIn contrast to the streamline approach, global tractography reconstructsall fiber tracts simultaneously, which enables fiber trajectories to be jointlyconsidered in determining the most possible fiber configuration. In this ap-proach, short fiber segments are bridged together to form the set of fibertracts that best explains the measured dMRI data. This globally optimaltract set is found by minimizing an energy function consisting of an internaland an external energy term: E(m) = Eint(m) + Eext(m,D). Here, m isthe set of all short fiber segments and D is the observed dMRI signal. Eintencourages the short fiber segments having consistent orientations to formchains and Eext forces the set of fiber tracts to be in agreement with D.We note that random changes are introduced into the reconstructed fiberset to avoid getting trapped in local minima. These random changes are ac-cepted or rejected through a Metropolis-Hastings sampler [42]. This globalapproach is particularly beneficial in resolving cases where local diffusioncharacteristics are ambiguous due to crossing fibers. Streamline tractogra-phy tends to terminate tracking at such locations, whereas global tractog-raphy exploits the geometry of surrounding fibers by considering the whole-brain fiber configuration in aggregate to compensate for the lack of reliablelocal information.2.2.3 ExperimentsMaterialsMinimally preprocessed RS fMRI and dMRI data from 40 unrelated healthysubjects of the Human Connectome Project (HCP) dataset [123] were used.2 The dataset comprised four RS fMRI scans of 15 minutes per subject,with a TR of 0.72 s and a voxel size of 2 mm (isotropic). The dMRI datahad a voxel size of 1.25 mm (isotropic), 3 shells (b=1000, 2000 and 3000s/mm2) and 288 gradient directions. Further details on acquisition can befound in [123].In addition to the minimal preprocessing that was already performed onthe HCP RS fMRI data which included gradient distortion correction, mo-tion correction, spatial normalization and intensity normalization [124], wefurther regressed out motion artifacts, white matter and cerebrospinal fluidconfounds, and principal components of high variance voxels found usingCompCor [125]. Band-pass filtering (0.01 Hz - 0.1 Hz) was subsequently ap-2This dataset is available online at: Estimating Anatomical Connectivityplied. To define ROIs, we functionally divided the brain into d=200 parcelsby temporally concatenating RS fMRI voxel time series across scans of allsubjects and applying Ward clustering [119]. This value of d provides higherfunctional homogeneity than typical anatomical atlases, while maintainingreasonable ROIs to time samples ratio for reliable connectivity estimation.The reason why we used a higher number of ROIs than that in Section2.1.4 is that HCP dataset has higher resolution, hence facilitates a morefine-grained parcellation. ROI time series were then generated by averagingvoxel time series within each ROI.The HCP dMRI data, which have been corrected for EPI distortion, eddycurrent, gradient nonlinearity and motion artifacts [124], were downsampledto 2 mm isotropic resolution to reduce the computational cost. Streamline(on both diffusion tensors and ODFs) and global tractography on ODFswere carried out using Dipy [122] and MITK [126], respectively. For allmethods, approximately 65,000 fibers per subject were reconstructed. Thefunctionally derived group parcellation map was warped onto the dMRIspace of each subject to facilitate the computation of fiber count.ResultsWe first present a qualitative comparison between streamline and globaltractography for a representative subject. Fig. 2.4 shows the tracts runningthrough an ROI within the corpus callosum as estimated using streamlinetractography on diffusion tensors and ODFs, and global tractography onODFs (in yellow, green and red, respectively). We chose this ROI since mosttractography algorithms have difficulties tracking the callosal projections[42]. As apparent in Fig. 2.4, streamline tractography only captured thedominant U-shaped callosal radiation. This result is expected for the caseof diffusion tensors since the principal diffusion directions of the tensors areambiguous at fiber crossings. However, the same pattern was observed evenwith diffusion modelled using CSA-ODF, which is able to delineate multipleintravoxel diffusion directions. This shows that streamline tractography isunable to make full use of the information in the ODFs. In contrast, globaltractography was able to better exploit the information in CSA-ODF andtracked the lateral transcallosal fibers.Quantitatively, the AC-FC correlation, computed by finding the Pear-sons correlation between intra-subject AC and FC estimates and averagingthese values across the group, was 0.2234±0.018 with streamline tractogra-phy on ODFs and 0.1643±0.012 with streamline tractography on diffusiontensors. The difference was found to be statistically significant at p<0.01362.2. Estimating Anatomical ConnectivityFigure 2.4: Fiber tracts reconstructed for a representative subject usingstreamline tractography on diffusion tensors (left), streamline tractographyon ODFs (middle), and global tractography on ODFs (right) with an ROIplaced in the corpus callosum (shown with a square). While streamline trac-tography results were limited to the callosal radiation, global tractographyalso tracked bilateral projections.based on the Wilcoxon signed rank test, showing that using ODFs can no-tably increase AC-FC consistency over diffusion tensors. With global trac-tography on ODFs, AC-FC consistency further increased to 0.2619±0.018.The difference between global and streamline tractography on ODFs wasstatistically significant at p<0.01. These findings confirm our hypothesisthat global tractography can better estimate the macroscale neural circuitry,leading to higher agreement between estimates of AC and FC, though thegain compared to using enhanced diffusion modelling, i.e. ODFs over ten-sors, was smaller. We speculate that uncertainty in fiber endpoint locationsis the limiting factor. Global tractography is driven by the measurements inareas with anisotropic diffusion, whereas diffusion around white-grey matterboundaries is near-isotropic. Thus, reconstructed fibers might not terminatein the correct location. We hence argue that modelling fiber endpoint uncer-tainty is an under-investigated problem that has high potential in improvingAC estimation.2.2.4 DiscussionStreamline tractography, which remains to be the most commonly used trac-tography approach, is known for its inability to resolve fiber crossings andpremature termination of the tracking process. On high quality multimodaldata, we demonstrated that this can lead to underestimation of the corre-372.3. Quantifying Anatomical Connectivitylation between estimates of AC and FC. We illustrated that this problemis more pronounced when streamline tractography is performed on conven-tional diffusion tensors as opposed to ODFs. We also showed that AC-FCconsistency can be improved by reconstructing the fiber tracts using globaltractography, which provides higher stability against noise and imaging ar-tifacts by jointly incorporating the geometry of all fiber tracts.Even though we showed that the consistency between AC and FC mea-sures can be improved by the choice of AC estimation method, higher cor-relation values have been previously reported in the literature [50, 63, 64].There are four possible reasons of this discrepancy. Firstly, we report ourvalues over all brain region pairs, instead of excluding absent anatomicalconnections as was the case in [63] or dividing the connection structure intosubsets as was done in [64]. Since anatomically unconnected brain regionpairs can still show high FC [97], analyzing only anatomically connectedbrain regions implies discarding part of AC-FC inconsistency, which can po-tentially inflate the AC-FC correlation values. Secondly, we do not applyany post-processing steps on the AC estimates and directly use fiber countsto assess the consistency between AC and FC measures. Fiber counts wereresampled into a Gaussian distribution in [63], which increased the averageAC-FC correlation from 0.156 to 0.36. Even though this strategy increasedthe consistency empirically, such resampling has no theoretical or physiolog-ical grounding. Thirdly, we do not apply distance correction to weigh downthe effect of longer anatomical connections in generating AC estimates aswas done in [50]. Such distance correction would lead to a positive bias inAC-FC correlation since proximal brain region pairs tend to have higher FCcompared to distal region pairs on average [97]. Finally, we do not applyglobal RS fMRI signal regression unlike [50, 63], since global signal regres-sion artificially introduces negative RS correlations [127]. We observed thatthe AC-FC correlation increased by more than 0.05 for both streamline andglobal tractography when the global signal was removed. However, whetherthis increase has a neural basis remains to be a point of debate.2.3 Quantifying Anatomical ConnectivitySection 2.3 is adapted from: M. Bajammal, B. Yoldemir, R. Abugharbieh. Compari-son of structural connectivity metrics for multimodal brain image analysis. In Proceedingsof the IEEE International Symposium on Biomedical Imaging (ISBI’15), pp. 934-937,Brooklyn, NY, USA, April 2015.382.3. Quantifying Anatomical ConnectivityMultimodal analyses typically require the quantification of connectivity, andinaccurate connectivity metrics may lead to flawed inferences. In contrastto quantification of FC, no well-analyzed and accepted metric exists forquantifying AC. To this end, we present a comparison of the commonly usedAC metrics assessed from the perspective of the largely accepted inherentrelationship between brain anatomy and function.As described in Section 1.4.1, the most commonly used AC metrics in-clude (i) fiber count, (ii) fiber count normalized by the total volume of theregion pairs they connect, (iii) total length of reconstructed fiber tracts,and (iv) average FA along fiber tracts connecting regions. Among these, ar-guably the most popular metric has been fiber count to date. We reiteratethat all of these AC metrics are confounded by several factors such as thenumber of seeds used for tracking the fibers and the structural features ofthe pathway (e.g. curvature, length and width). Nonetheless, these metricsare increasingly used in the literature driven by a practical need for quan-tifying the degree of connection between brain regions. We thus perform aquantitative analysis to determine which AC metric has the highest poten-tial of being of practical use in multimodal brain image analysis efforts, withthe understanding that these AC metrics do not truly convey the “strength”of the anatomical connections.2.3.1 Anatomical Connectivity MetricsLet rki,j be the kth reconstructed fiber between a pair of anatomically con-nected regions Pi and Pj . We consider four widely used anatomical connec-tivity measures in this thesis: fiber count (fi,j), fiber count normalized bythe total volume of the connected regions (Ni,j), total length of fibers con-necting region pairs (Li,j), and average FA along the fibers. For each subject,we compute the Pearson’s correlation between the FC and AC estimates toquantify the AC-FC relationship for each AC metric. More formally, themetrics can be expressed asNi,j =fi,jV (Pi) + V (Pj)(2.6)Li,j =∑k`(rki,j), (2.7)where V (·) is the volume of the corresponding region, and `(rki,j) is the lengthof rki,j .Prior to the computation of AC metrics, we reconstruct the fibers viaglobal tractography on CSA-ODF using MITK [126]. Global tractography392.3. Quantifying Anatomical Connectivitywas chosen over the more common streamline tractography since we haveshown in Section 2.2 that it facilitates higher AC-FC consistency.2.3.2 ExperimentsMaterialsdMRI, RS fMRI and task fMRI data from HCP dataset were used in the ex-periments [123]. Details and preprocessing steps of dMRI and RS fMRI dataand the generation of the brain parcellation map are presented in Section2.2.3. Each subject had 7 task fMRI datasets (two 2-5 minute scans eachfor working memory, gambling, motor, language, social cognition, relationalprocessing, and emotional processing) obtained using the same acquisitionparameters as the RS fMRI data. Preprocessing of the task fMRI dataclosely followed that of RS fMRI, except a high-pass filter at 1/128 Hz wasused instead of a bandpass filter.ResultsThe results of the AC-FC correlation for RS and task fMRI are shown inFigs. 2.5 and 2.6. As observed from these figures, average FA has lowercorrelation with FC compared to the rest of the examined AC metrics (fibercount, volume-normalized fiber count, and total fiber length). This is truefor both RS and task fMRI. We speculate that the observed low average FAcorrelation can be attributed to the large number of factors affecting localdiffusion anisotropy [53]. These results also show that volume-normalizedfiber count has the highest correlation with FC compared to the rest of theexamined AC metrics for both RS and task fMRI. The pairwise differencesbetween AC-FC correlation assessed using normalized fiber count and otherAC metrics were found to be statistically significant at p<0.001 based on theWilcoxon signed rank test. Our results thus imply that the compensationfor the differences in number of fibers due to the variable size of brain regionsyields better depiction of anatomical networks.We also note that the relatively consistent AC-FC correlation levelsacross a variety of tasks and RS data support the notion that AC formsthe backbone of the brain connectivity around which functional reorganiza-tion occurs to respond to different tasks. A diverse repertoire of functionalbrain connectivity patterns can thus arise constrained by the same anatom-ical substrate.AC-FC correlation is not the only factor to consider while choosing anAC metric. It is known that proximal brain region pairs tend to have higher402.3. Quantifying Anatomical Connectivity5 10 15 20 25 30 350.  fiber countnormalized fiber countaverage FAtotal fiber lengthAC−FC CorrelationSubjectFigure 2.5: Correlation between AC and FC using four commonly usedAC metrics. FC was estimated from RS fMRI data, thus depicts intrinsicconnections.5 10 15 20 25 30 3500.  fiber countnormalized fiber countaverage FAtotal fiber lengthAC−FC CorrelationSubjectFigure 2.6: Correlation between AC and FC using four commonly used ACmetrics. FC was estimated from task fMRI data acquired during 7 differenttasks. The narrow shaded bands in represent the standard deviation ofAC-FC correlation across different tasks.412.3. Quantifying Anatomical ConnectivityTable 2.1: Comparison of fiber length bias for AC metrics.AC metric Fiber length biasFiber count 0.14Normalized fiber count 0.17Average FA 0.15Total fiber length 0.10FC compared to distal region pairs on average [97]. An AC metric consider-ing mostly shorter fibers and ignoring the long fibers might thus lead to aninflated AC-FC consistency, even though there is no reason to believe longfibers should be weighed less than short fibers when quantifying AC. Indeed,AC metrics have an inherent bias in favor of shorter fibers since shorter fibersare easier to be reconstructed. This in turn induces a dependency of the AC-FC relationship on distance. To quantify the extent of this phenomenon, wecomputed the difference between group-averaged AC-FC correlation valuesof the longest and shortest 20% of all fibers, which we refer to as fiber lengthbias. As shown in Table 2.1, total fiber length metric has the smallest fiberlength bias. This finding confirms that total fiber length indeed partiallycorrects for length bias fiber count metrics suffer from. However, this metricresults in a considerably lower AC-FC correlation than fiber count and nor-malized fiber count. We thus used fiber count to quantify AC in this thesis,which provides a lower fiber length bias than normalized fiber count, and ahigher AC-FC correlation than total fiber length.2.3.3 DiscussionWe presented a comparison of four AC metrics computed from tractogra-phy results with respect to their relationship to FC. Among the metricsconsidered, we showed that volume-normalized fiber count has the high-est correlation with FC for both RS and task data. On the other hand,our results showed that average FA has the lowest correlation with FC. Wespeculate that the reason of this low correlation is the non-specificity of FA,with several inadvertent factors (such as axonal density, axonal ordering,and amount of myelination) modulating it along with the reliability of localdiffusion information [53]. In addition, we also demonstrated that total fiberlength metric reduces the fiber length bias associated with shorter fibers. Tostrike a balance between high AC-FC consistency and low fiber length bias,we chose to employ fiber count to quantify AC in this thesis.422.4. Summary2.4 SummaryIn this chapter, we investigated several steps leading to quantitative AC es-timates, particularly from the perspective of their effect on the consistencybetween AC and FC. We showed that part of the inconsistency between ACand FC estimates can be attributed to the analysis methods used to obtainmeasures of AC. In particular, we showed that increasing the spatial res-olution of DW images using a dictionary based super-resolution approachimproves the reconstructed fiber tracts and tract-density maps qualitatively,and leads to improved AC-FC consistency. Further, we showed that mod-elling local diffusion using ODFs over diffusion tensors, and reconstructingthe fibers using global tractography over streamline tractography help betterresolve fiber crossings around corpus callosum, and provide higher AC-FCconsistency over the whole brain. Finally, we showed that quantifying ACusing fiber count or volume-normalized fiber count results in higher AC-FCconsistency compared to average FA and total fiber length. Our experi-ments revealed that fiber count has less bias in favor of shorter fibers thanvolume-normalized fiber count, and thus fiber count is the metric used toquantify AC in this thesis.43Chapter 3Multimodal Brain ActivationDetectionTo map brain areas to function, the standard analysis approach modelsfMRI observations as a combination of expected temporal responses usingGLM [26]. However, the strong noise in fMRI data arising from confounds,such as scanner drifts, motion artifacts, and physiological effects, greatlyhampers reliable detection of brain activation. Exploring other sources ofextra information to regularize activation detection has thus proven to bebeneficial. To model the integrative property of the brain, which is knownto facilitate brain function [128], the use of local neighbourhood informationhas been proposed to regularize activation detection [69, 72]. Although suchmethods help suppress false spatially-isolated activations by encouragingspatial continuity, they completely ignore long-range functional interactions.The incorporation of functional connectivity information into task activationdetection has also been put forth [74], but the approach taken estimates bothactivation effects and functional connectivity from the same dataset, hencethe information gain might be limited. Other works investigated the use ofRS functional connectivity information to inform task activation detection[79] as motivated by the similarity of RS subnetworks and those engagedduring task [78]. Adding to the arsenal of prior-informed brain activationdetection methods, we propose two novel multimodal activation detectionapproaches in this chapter.Firstly, we explore the feasibility of deconfounding the effects of RS ac-tivity from task fMRI data. Recent studies suggest that part of the noise intask fMRI data actually pertains to ongoing RS brain activity. Due to thesporadic nature of RS temporal dynamics, pre-specifying temporal regressorsto reduce the confounding effects of RS activity on task activation detectionis far from trivial. We propose a novel approach that exploits the intrinsictask-rest relationships in brain activity for addressing this challenging prob-lem. With an approximate task activation pattern serving as a seed, we firstinfer areas in the brain that are intrinsically connected to this seed using RSfMRI data. We then apply PCA to extract the RS component within the443.1. Modelling Intrinsic Functional Connectivitytask fMRI time courses of the identified intrinsically-connected brain areas.Using the learned RS modulations as confound regressors, we re-estimatethe task activation pattern, and repeat this process until convergence. Onreal data, we show that removal of the estimated RS modulations from taskfMRI data significantly improves activation detection (Section 3.1).Secondly, we propose incorporating anatomical connectivity into brainactivation detection as motivated by how the functional integration of dis-tinct brain areas is facilitated via neural fiber pathways. We formulate acti-vation detection as a probabilistic graph-based segmentation problem withfiber networks estimated from dMRI data serving as a prior. Our approachis reinforced with a data-driven scheme for refining the connectivity priorto reflect the fact that not all fibers are necessarily deployed during a givencognitive task as well as to account for false fiber tracts arising from limita-tions of dMRI tractography. Validating on real clinical data collected fromschizophrenia patients and matched healthy controls, we show that incorpo-rating anatomical connectivity significantly increases sensitivity in detectingtask activation in controls compared to existing univariate techniques. Fur-ther, we illustrate how our model enables the detection of significant groupactivation differences between controls and patients that are missed withstandard methods (Section 3.2).3.1 Modelling Intrinsic Functional ConnectivityThe frequency range at which RS activity resides is typically found to bebetween 0.01 to 0.1 Hz [76], which overlaps with the stimulus frequenciesemployed in most task fMRI studies. Thus, standard high pass filtering 1/128 Hz, which is the default cutoff frequency in the SPM software, forremoving temporal drifts in task fMRI data, would not account for ongoingRS modulations. Also, unlike task-evoked responses, which are time-lockedto stimulus, the time at which RS activity peaks and troughs is difficult topredict. Hence, pre-specifying temporal regressors to model RS activity isnon-trivial.Albeit its seemingly sporadic temporal dynamics, RS activity is not ran-dom [76]. Rather, strong synchrony in RS modulations between specificbrain areas has been observed in numerous studies [76]. In fact, many ofSection 3.1 is adapted from: B. Yoldemir, B. Ng, R. Abugharbieh. Deconfoundingthe effects of resting state activity on task activation detection in fMRI. In Proceedingsof the MICCAI Workshop on Multimodal Brain Image Analysis, LNCS, vol. 7509, pp.51-60, Nice, France, October 2012.453.1. Modelling Intrinsic Functional Connectivitythe detected RS subnetworks exhibit high resemblance to subnetworks seenin task experiments [78]. Further supporting this finding is a recent work[79] that demonstrated enhanced sensitivity in task activation detection byincorporating an RS connectivity prior. In addition, studies that jointly ex-amined RS fMRI and dMRI data indicate an anatomical basis for RS activity[97, 98]. Thus, the spatial patterns of RS subnetworks would presumablybe constrained by the underlying fiber pathways [97, 98]. Taken together,these findings suggest that there is spatial structure in RS activity and thatthe spatial structure of ongoing RS activity during task would likely remainsimilar to that during rest [75].In this section, we propose a novel approach for RS activity removal ina general setting. The key challenge to this problem is that RS activity isinternally-driven by the brain, as opposed to being evoked by external stimu-lus with known timing. It is thus not obvious how the temporal dynamics ofRS modulations that occurred during task performance can be determineda priori. Representative time courses reflective of ongoing RS activity musthence be extracted from the task fMRI data itself. Since the brain com-prises multiple subnetworks [78], the RS modulations superimposed on thefMRI responses of the task-activated brain areas would be specific to the RSsubnetwork in which these brain areas belong. Extracting RS activity fromtask fMRI data would thus require knowing the parts of the brain that areactivated and their intrinsically-connected areas, which introduces a circularproblem. To deal with this issue, we employ an iterative strategy in whichwe first apply seed-based analysis [22] with an approximate task activationpattern being the seed to infer the intrinsically-connected brain areas fromRS fMRI data. Assuming the spatial structure of RS subnetworks is sus-tained during task performance [75], we extract RS modulations from thetask fMRI time courses of the identified brain areas and re-estimate the taskactivation pattern with the learned RS activity as confound regressors. Theoverall method is summarized in Fig. Group Activation DetectionOur approach begins with the estimation of an approximate task activationpattern (i.e. without accounting for ongoing RS modulations). A standardGLM is first applied to compute the intra-subject activation effects [26]:Yi = Xiβi + Ei, (3.1)where Yi is a t× d matrix containing the task fMRI time courses of d brainareas of subject i, Xi = [Xtask|Xiconfounds] is a t × p matrix with Xtask463.1. Modelling Intrinsic Functional ConnectivityGroup activation detectionRS subnetwork detectionRS activityestimationFigure 3.1: Graphical depiction of proposed RS activity removal approach.Xi = [Xtask|Xiconfounds] is a regressor matrix, where Xtask corresponds totask regressors and Xiconfounds corresponds to confound regressors specificto subject i. Yi are the task fMRI time courses of subject i. ΛA is the setof activated brain areas common across a group of subjects. Zi are the RSfMRI time courses of subject i. ΛiC is the set of brain areas estimated to beintrinsically connected to ΛA for subject i. yiRS is the estimated RS activitytime course of subject i, which is entered into Xi as a confound regressorfor re-estimating the group activation pattern ΛA. The three steps: groupactivation detection, RS subnetwork detection, and RS activity estimation,are repeated until ΛA stabilizes.corresponding to task regressors and Xiconfounds corresponding to confoundregressors specific to subject i, βi is a p × d activation effect matrix tobe estimated, and Ei is a t × d residual matrix. Due to the strong noisein task fMRI data, activation patterns estimated at the intra-subject levelmight be inaccurate [129]. Therefore, we opt to combine information acrosssubjects in generating a group activation map, which is then used as a seedfor identifying intrinsically-connected brain areas. To infer group activation,we apply a max-t permutation test [67] on βi of all subjects, which implicitlyaccounts for multiple comparisons and provides strong control over falsedetections. Group activation is declared at a p-value threshold of 0.05. Wedenote the set of detected brain areas as ΛA.3.1.2 RS Subnetwork DetectionWith the detected group activation pattern taken as a seed, our goal is toidentify brain areas that belong to the same RS subnetwork as the seed,so that we can extract RS modulations specific to this RS subnetwork. Toproceed, we first average the RS fMRI time courses within the detectedactivated brain areas in generating a seed time course for each subject.473.1. Modelling Intrinsic Functional ConnectivityWe then apply the standard seed-based analysis [22] to find brain areasintrinsically-connected to this seed for each subject i:Zi∼S = ziswi + Ωi, (3.2)where Zi∼S is an n×(d−|ΛA|) matrix containing the RS fMRI time courses ofall brain areas except those in ΛA, |ΛA| is the number of brain areas in ΛA,zis is an n× 1 vector containing the seed time course, wi is a 1× (d− |ΛA|)vector with each element reflecting the correlation between the seed andeach brain area not in ΛA, and Ωi is an n× (d−|ΛA|) residual matrix. Sta-tistical significance of each element of wi is declared at a p-value thresholdof 0.05 with false discovery rate (FDR) correction [130] to account for multi-ple comparisons. FDR correction is used instead of max-t permutation testdue to correlations between brain volumes at adjacent time points, whichviolates the independent sample assumption in max-t permutation test [67].A max-t permutation test could be applied to identify brain areas that aresignificantly correlated with the seed at the group level. However, comparedto task-based experiments, RS experiments are less prone to motion arti-facts, which constitute a major part of fMRI noise. Reliable RS subnetworkscould thus potentially be extracted at the intra-subject level. We hence optto perform intra-subject seed-based analysis to retain subject-specific infor-mation. We denote the set of brain areas significantly correlated with theseed as ΛiC .3.1.3 RS Activity Estimation and RemovalAfter finding the set of brain areas that is intrinsically-connected to theestimated task activation pattern for each subject i, the next step is toextract the RS components from the task fMRI time courses of these brainareas, which we denote as YiC . To target the specific frequency range atwhich RS activity resides, we first band-pass filter each column of YiC atcut-off frequencies of 0.01 and 0.1 Hz. Since the estimated task activationpattern is only an approximation without accounting for the confoundingeffects of RS activity, some intrinsically-connected brain areas might in factbe activated considering the resemblance between task and RS subnetworks[78]. Thus, YiC might contain task signals. To remove the task-relatedresponse in YiC , we apply PCA through eigen-decomposition to separateYiC into task and non-task components:Ci = UiDiUiT , (3.3)483.1. Modelling Intrinsic Functional Connectivitywhere Ci is the d×d covariance matrix of YiC , Ui is a d×d matrix containingthe eigenvectors of Ci, and Di is a d× d matrix containing the eigenvaluesof Ci along the diagonal. The columns of Ui are ordered such that the firstcolumn, Ui1, corresponds to the largest eigenvalue. To identify the task-related components, we compute the correlation between each column of Uiand the task regressor. Statistical significance in correlation is declared ata p-value threshold of 0.05 with FDR correction. For the data examinedin this work, the task regressor is found to be most significantly correlatedwith Ui1. This high correlation between Ui1 and the task stimulus, as shownin Fig. 3.2, signifies a definite need for task response removal from YiC . Weremove the task components by reconstructing YiC with the significantlycorrelated columns of Ui discarded.  Task regressorMean first PCAmplitudeTime indexFigure 3.2: Dominant principal component (PC) extracted from time coursesof intrinsically-connected brain areas. The thick red and blue lines corre-spond to task regressor and mean dominant PC across subjects. Each thinline corresponds to the dominant PC of a single subject.Denoting the reconstructed YiC as Vi∼task, we take the mean of Vi∼taskover brain areas to generate a representative RS activity time course foreach subject i, yiRS , which we enter into (3.1) as a confound regressor tore-estimate the task activation pattern. This process is repeated until thegroup activation pattern ΛA stabilizes.493.1. Modelling Intrinsic Functional Connectivity3.1.4 ExperimentsMaterialsFor testing our proposed RS activity removal approach, we used the publiclyavailable Multiband Test-Retest Pilot Dataset, which was released as a partof the 1000 Functional Connectomes Project. 3 Excluding subjects withmissing brain volumes, the dataset comprises 19 subjects (14 men, 5 women,mean age 33.1±13.2 years). Each subject performed a passive viewing taskin which a checkerboard was displayed on a monitor for 20 s, with 20 sof rest interleaved between stimulus blocks. The total task duration wasapproximately 2.5 minutes. Task fMRI data were acquired with a TR of 1.4s and a voxel size of 2 mm (isotropic). RS fMRI data of 5 minutes durationwere also collected with a TR of 2.5 s and a voxel size of 3mm (isotropic).Preprocessing of the data followed the pipeline presented in Section 2.2.3.ResultsTo validate our proposed approach, we compared applying the standardunivariate analysis [26] with and without RS activity removal in detectinggroup activation. We denote these two cases as RSR and GLM, respectively.For increased group activation detection to be a legitimate validation cri-terion, strong control over false positive rate is critical. For this, we usedthe max-t permutation test [67] for both RSR and GLM, which implicitlyaccounts for multiple comparisons, provides strong control on false posi-tive rate, and generates less conservative t-value thresholds than Bonferronicorrection [67].Fig. 3.3 shows the number of detected parcels for different p-valuethresholds. For the same specificity, our approach provided higher detectionsensitivity than GLM in general. To assess whether the increased detec-tion was statistically significant, we employed a “parcel-label” permutationtest. Specifically, for each permutation, we first randomly selected half ofthe parcels and exchanged the labels (i.e. active or non-active) assigned byRSR and GLM for each p-value threshold. We then computed the differ-ence in the number of detected parcels with and without RS removal, whichwe denote as Ndiff . This procedure was repeated 1000 times to generatea null distribution. The original Ndiff was found to be greater than the95th percentile of the null distribution for all corrected p-value thresholdswithin the typical range of [0.01, 0.05]. The detection improvement with3This dataset is available online at: Modelling Intrinsic Functional ConnectivityRSR compared to GLM was thus statistically significant. We note that im-provement in detection was observed even with just one iteration of RSR,and the detected activation pattern stabilized within two iterations, i.e. nomore than a couple of parcels changing labels in subsequent iterations.0 0.01 0.02 0.03 0.04 0.0505101520253035GLMRS Removed-Iter#1RS Removed-Iter#2RS Removed-Spat. Rand.RS Removed-Temp. Rand.p-value threshold (corrected)# of activated parcelsFigure 3.3: Comparison of activation detection results with and without theproposed RS activity removal method. Number of parcels detected withsignificant activation vs. p-value thresholds.To show that there is specific temporal structure in the estimated RSactivity time courses that gave rise to the observed detection improve-ment, we applied RSR on temporally-permuted RS activity time courses100 times. The average number of detected parcels over the 100 permuta-tions (RS Removed-Temp. Rand. in Fig. 3.3) was found to be similar tothat of GLM. The difference in detection performance between RSR andtemporally-permuted case was statistically significant based on the parcel-label permutation test with a threshold set at the 95th percentile of the nulldistribution. Our results thus indicate that there is certain temporal struc-ture in the estimated RS activity time courses that is critical for successfulRS activity removal.Since the brain comprises multiple subnetworks [78], not all parcels wouldcontain the same RS modulations as those superimposed on the underlyingtask activated brain areas. To illustrate this point, we applied RSR withRS modulations extracted from Nc randomly selected parcels (excludingparcels identified by RSR), where Nc is the number of intrinsically-connected513.1. Modelling Intrinsic Functional Connectivitybrain areas originally determined with RSR. The average number of detectedparcels over 100 random subsets of parcels (RS Removed-Spat. Rand. in Fig.3.3) was found to be similar to that of GLM for p-value thresholds between 0and 0.02, and modestly better than GLM for p-value thresholds above 0.02.We suspect the increased detection arises from how some parcels might beintrinsically-connected to the task-evoked brain areas, but the estimatedcorrelations were declared not significant due to noise. Such parcels wouldcontain RS modulations common to the task activation pattern, hence theincreased activation detection observed. Nevertheless, the increase was notstatistically significant based on the parcel-label permutation test with athreshold set at the 95th percentile of the null distribution.Qualitatively, RSR additionally detected brain areas adjacent to thosefound by GLM (Fig. 3.4). More bilateral activation was also found withRSR. The detected brain areas lie within the primary visual cortex andthe extrastriate cortex, which are known to pertain to visual checkerboardstimulus [131], hence confirming our results.Figure 3.4: Detected activation patterns with and without the proposed RSactivity removal method. Three axial slices shown. Parcels detected at p-value < 0.05 (corrected). Red = detected by RSR only. Purple = detectedby both GLM and RSR.3.1.5 DiscussionWe proposed a novel approach for the estimation and removal of continualRS activity in task fMRI data. Exploiting how the spatial structure ofRS subnetworks are constrained by the underlying fiber pathways hencewould remain similar during task performance, our approach first extractsRS modulations from task fMRI time courses within brain areas that are523.2. Modelling Anatomical Connectivitysignificantly correlated with an approximate task activation pattern at rest.The estimated RS modulations are then entered into a GLM as confoundregressors to model the effects of RS activity. Applying our approach onreal data resulted in statistically significant improvement in task activationdetection. In agreement with the seminal work by Fox et al. [75], ourresults indicate that RS activity also contributes to the noise seen in taskfMRI data, in contrast to the traditional belief that only scanner artifacts,head motions, and physiological confounds contribute to fMRI noise. It isthus important to model RS activity in task fMRI studies.3.2 Modelling Anatomical ConnectivityIn this section, we propose incorporating AC information into task activa-tion detection. Given that fiber pathways serve as the physical substrate forfunctional interactions, we hypothesize that intrinsically connected brain ar-eas would likely be in similar state, e.g. co-activate, during task [79]. Thus,informing activation detection with AC should presumably improve the de-tection sensitivity. For this, we employ the graph-theoretic random walker(RW) formulation [132], which easily permits such an integrated scheme forestimating activation probabilities. Posterior activation probabilities esti-mated by the RW formulation are guaranteed to be unique and globally-optimal [132], which makes RW an eminent choice. RW has been previouslyapplied to task activation detection with functional connectivity taken asthe prior [74]. Here, we investigate the implications of complementing taskactivation detection analysis with AC information. To infer group activationfrom posterior activation probabilities, we devise a permutation test withactivation probabilities as attributes, which we empirically show to providestronger control on false positive rate than simply comparing the posteriorprobability of being activated, not activated, and de-activated. After es-timating a group activation pattern, we iteratively refine the anatomicalconnectivity prior by removing the links between non-active brain regionsand re-estimate the posterior activation probabilities until the detected ac-tivation pattern stabilizes. An overview of our multimodal task activationdetection approach is shown in Fig. 3.5.On real data, we demonstrate that incorporating AC increases sensi-Section 3.2 is adapted from: B. Yoldemir, B. Ng, T. Woodward, R. Abugharbieh.Fiber connectivity integrated brain activation detection. In Proceedings of the 23rd Bien-nial International Conference on Information Processing in Medical Imaging (IPMI’13),LNCS, vol. 7917, pp. 135-146, Asilomar, CA, June 2013.533.2. Modelling Anatomical ConnectivityAnatomical connectivity prior estimationActivation likelihood estimationPosteriorprobability estimationGroupactivationinferencedMRI datafMRI dataActivation statesFigure 3.5: Overview of the proposed multimodal task activation detec-tion approach. Anatomical connectivity prior and activation likelihood esti-mates, i.e. label priors, are integrated under the RW formulation to find theposterior activation probabilities. A permutation test is applied on theseprobabilities to infer group activation. The anatomical connectivity prior isthen iteratively refined based on the activation states of the brain regionsuntil convergence.tivity in detecting group activation over using univariate techniques. Wefurther show that our method is able to detect significant group activationdifferences between schizophrenia patients and healthy controls, which aremissed by standard analysis approaches.3.2.1 Random Walker for Activation EstimationRW is a graph-based image segmentation approach, in which graph nodes(vertices) correspond to image voxels and graph edges connecting neigh-bouring nodes are assigned weights reflecting node similarity. In its originalformulation [133], RW labels nodes based on the probability that a randomwalker starting from each node will first reach a pre-labeled seed given edgeweights that bias the paths. This requires user-specified seeds and does notutilize local observations at each vertex. Thus, in the context of activationdetection with brain regions being the graph nodes and their interactionsmodelled through graph edges, this would require user-defined seeds for ev-ery functionally-disparate region and only functional interactions would beconsidered without accounting for activation effects. We thus adopt theformulation in [132], which overcomes the need for user interaction and in-tegrates activation effects into the formulation as label priors. This is equiv-alent to adding a floating vertex for each label and connecting these floatingvertices to every vertex in the original graph with label priors being the543.2. Modelling Anatomical Connectivityweights of the added edges [132]. In this formulation, posterior probabilitiesare calculated by the minimization of the following energy functional:E(xs) = xsTLxs +K∑k=1,k 6=sxkTΛkxk + (xs − 1)TΛs(xs − 1), (3.4)where xs is an M × 1 vector of unknown posterior probabilities of each ROIbelonging to class s, L is an M×M weighted Laplacian matrix, Λs is an M×M diagonal matrix having the prior probabilities of the ROIs belonging toclass s on its diagonal, K is the number of class labels, and M is the numberof brain regions. The first term in (3.4) is a spatial term for modelling theinteractions between graph vertices as characterized by L. The second termdenotes the aspatial component for modelling the local observations at thevertices. In our context of activation detection, the spatial term modelsthe anatomical connectivity information and the aspatial term models theactivation effects. The method can thus be thought of as grouping brainregions into classes via random walk on an augmented graph, where edgesin the original graph are weighted by AC information and edges leadingto the floating nodes are weighted by activation effects. Assuming equalweighting between the spatial and aspatial energy terms in (3.4), it hasbeen shown [132] that the posterior probabilities can be found by solving(L +K∑k=1Λk)xs = λs, (3.5)where λs is an M×1 vector consisting of the diagonal elements of Λs. SinceL is positive semi-definite and Λk is strictly positive definite, their summa-tion would be diagonally dominant. Hence, matrix inversion is possible forsolving (3.5). Following [134], we set K = 3 and define the class labels as de-active, nonactive, and active. For clarity, we explicitly denote the posteriorprobabilities for each class as pD, pN and pA, corresponding to deactive,nonactive and active classes, respectively.3.2.2 Anatomical Connectivity Prior EstimationLet D be an M ×M weighted adjacency matrix, where each element Dij isan estimate of the AC between brain regions i and j, set to fiber count fol-lowing our discussion in Section 2.3. The corresponding M ×M normalized553.2. Modelling Anatomical ConnectivityLaplacian matrix, L, of D is given byLij =1 if i = j and dj 6= 0− 1√didjif i and j are adjacent0 otherwise, (3.6)where di =∑j Dij is the degree of node i. A major difficulty with tractog-raphy is resolving fiber crossing regions where accuracy of most algorithmsis seriously affected. In [64], it was proposed that multiplying D by itselfmay help address this problem by generating multi-step fibers from parts offibers that might be split at crossing regions:DMS = exp(D) =∞∑k=01k!Dk, (3.7)where exp(·) denotes the matrix exponential. Dk = D ·D ·D · ... and Dkijis the number of paths of length k connecting regions i and j. DMS hencecomprises all possible paths between each region pair, where indirect pathsare more heavily penalized as these paths are potentially artifactual [64].3.2.3 Activation Likelihood EstimationTo estimate the activation likelihoods, which are used as label priors λsin RW, we first apply the classical GLM [10] to compute the intra-subjectactivation statistics:yj = Xβj + ejtj = βˆj/se(βˆj),(3.8)where yj is the n × 1 time course of ROI j, X is the n × r design matrixof expected responses, βˆj is an r × 1 vector containing estimates of theactivation effects, βj , ej is the n×1 residual assumed to be white Gaussiannoise after preprocessing, se(βˆj) is the standard error of βˆj , tj is the r × 1vector of sought activation statistics, n is the number of time points, and ris the number of experimental conditions. Columns of the design matrix Xare generated by convolving the canonical hemodynamic response function(HRF) with a boxcar time-locked to stimulus [10]. To compute the priorprobabilities of ROIs belonging to each class, we fit a Gamma-Gaussian-Gamma (GGG) mixture model separately to the t-values of each condition,tc, i.e. cth element of tj assembled across all ROIs [134]. The Gaussian563.2. Modelling Anatomical Connectivitydistribution models the nonactive state and the Gamma distributions modelthe deactive and active states:tcj ∼ piDΓ(kD, θD) + piNN (µ, σ) + piAΓ(kA, θA), (3.9)where µ is the mean and σ is the standard deviation of the Gaussian com-ponent, kD and kA are the shape parameters, and θD and θA are the scaleparameters of the deactivation (D) and activation (A) components, respec-tively. We employ the expectation-maximization (EM) algorithm [135] toestimate the model parameters separately for each experimental conditionc with the probabilities of tc given the parameter estimates used as the la-bel priors λs. These label priors and the Laplacian matrix are combinedthrough (3.5) to estimate the posterior activation probabilities.3.2.4 Group Activation InferenceThe high dimensionality of fMRI data elicits a high risk of false detection.To infer group activation from activation statistics, such as t-values, severalmethods that control for false positive rate have been proposed, e.g. Bon-ferroni correction, Gaussian random field theory, max-t permutation test[67]. For group activation inference from posterior activation probabilitymaps, most studies directly threshold the posterior activation probabilitiesat 1/K, where K is the number of classes, arguing that activation inferencefrom posterior activation probabilities does not suffer from false positives[136]. However, as we will empirically demonstrate in the next subsection,directly thresholding the posterior probabilities is actually prone to falsedetection, necessitating a more rigorous activation inference method. Tothis end, we propose a permutation test for group inference from activationprobabilities that controls for the false positive rate. Specifically, for eachpermutation, we first randomly select one third of the subjects and swap theposterior probabilities pA and pD of each selected subject. We note thatthis swap is done at the intra-subject level for all ROIs, hence the spatialpattern of the activation probabilities is preserved. Similarly, we swap pAand pN for another third of randomly selected subjects. We then computethe z-scores of pA for each permutation across all subjects. This procedureis repeated 10,000 times to generate the null distribution of activation prob-abilities of each ROI. We assign the p-value for the activation likelihood ofeach ROI as the number of times the z-scores of permuted pA are greaterthan the z-scores of the original pA values divided by the total number oftrials (in this case, 10,000). Under the null hypothesis that there is no ac-tivation, the z-score of the original pA value of an ROI would lie around573.2. Modelling Anatomical Connectivitythe mean of the generated null distribution, resulting in a non-significant p-value around 0.5. Finally, FDR [137] is applied to these p-values to accountfor multiple comparisons. As shown in the next subsection on a syntheticcase, this permutation test offers a much lower false positive rate comparedto posterior probability thresholding while keeping the true positive rate atthe same (or even at a higher) rate.3.2.5 Re-estimation of Group Activation Using PartialConnectomeSince some of the estimated fiber tracts might be false due to tractographyerrors and not all fibers are necessarily employed during a given cognitivetask, we propose restricting the anatomical connectivity prior to the sub-set of fibers that are likely to be active by iteratively refining the originalconnectivity prior, Dorig, asDm+1 = Dtask +1log(m+ 1)Dm (3.10)Dtaskij ={Dorigij if `mi ∨ `mj = 10 if `mi = `mj = 0, (3.11)where Dm is the M × M partial adjacency matrix found in iteration mwith D1 set to Dorig. `mi is the estimated activation state of ROI i duringiteration m with a value of 1 denoting activated and 0 denoting not activatedor deactivated. We note that updating Dorig using the above scheme enablesDorig to gradually evolve as opposed to having its information from previousiterations completely discarded. The overall process proceeds by estimatingthe group activation map with D1 as the prior. We then refine D1 using(3.10) and (3.11), and repeat the process until the group activation mapstabilizes.3.2.6 ExperimentsMaterialsAfter obtaining informed consent, fMRI data were collected from 13 healthysubjects (6 men, 7 women, mean age 27.46± 6.38 years) and 7 schizophreniapatients (5 men, 2 women, mean age 30.57± 10.08 years). Each subject wasfirst presented with words in four different contexts: associating, hearing,solving and reading, during a non-scanned encoding session. The subjects583.2. Modelling Anatomical Connectivitywere then presented the same set of words in a subsequent recall run duringwhich fMRI data were acquired and subjects were asked to indicate thecontext in which the presented words were previously encountered. Imageacquisition was performed on a Philips Achieva 3.0 T MRI scanner usinga T2*-weighted gradient-echo spin pulse sequence with a repetition time of2000 ms, an echo time of 30 ms, a flip angle of 90◦, a field of view of 240×240mm, and an in-plane resolution of 80×80 pixels. Each volume comprised 36axial slices of 3 mm thickness with a 1 mm gap. Each scan lasted for 920 s,which tallies to 460 fMRI volumes.dMRI data were collected from the same subjects using a Philips Achieva3.0 T MRI scanner with a TR of 7500 ms, a TE of 54 ms, an EPI factorof 59, an FOV of 224×224 mm and an in-plane resolution of 256×256 pix-els. Fifteen diffusion weighted volumes were acquired at a b-value of 800s/mm2 in addition to a volume with no diffusion sensitization. Each volumeconsisted of 72 slices of 2 mm thickness with no gap. Acquisition time was480 s. Preprocessing of the fMRI and dMRI data followed the pipeline pre-sented in Section 2.2.3 and Section 2.3.2. We used MedINRIA [138] for fibertractography. To facilitate the computation of fiber count, we warped ourfunctionally derived group parcel map to the b=0 volume of each subject.ResultsWe first present the synthetic test performed to assess our proposed permu-tation test as compared to the posterior probability thresholding approachcommonly employed in the literature [136]. For evaluation on real data, wecompare the sensitivity of our proposed approach in detecting group activa-tion in controls against that of univariate techniques. We then contrast ourmethod against classical schemes in detecting group activation differencesbetween schizophrenia patients and controls.We performed a synthetic test carefully designed from activation proba-bilities calculated from the real data of healthy controls. Specifically, afterestimating the group activation map with our approach, we used the high-est one third of pN among nonactive ROIs and the highest one third ofpA among active ROIs to generate a pseudo ground truth of nonactivationand activation probability distributions at the intra-subject level. Thesepseudo ground truth distributions are assumed to be Gaussian with meansand standard deviations set to that of the respective thirds of pN and pA.Out of a total of 100 synthetic datasets, each dataset comprised 13 subjectshaving 500 ROIs each, with 100 of them defined to be active. Randomsamples of pN and pA were drawn for each subject from the corresponding593.2. Modelling Anatomical Connectivityprobability distributions. Deactivation probabilities were computed basedon how posterior probabilities should sum to 1. Assessing group activationwith the proposed permutation test resulted in a true positive rate (TP) of0.819±0.038 and a false positive rate (FP) of 0.008±0.004, whereas thresh-olding pA at 1/3 gave a TP of 0.374±0.049 and an FP of 0.097±0.015.One-sample t-tests among TPs and FPs of these two strategies declared thedifferences to be significant at p< 10−6, demonstrating superior sensitivityand specificity.We compare the sensitivity of our method against classical GLM andagainst inferring group activation from the activation likelihoods given bythe GGG mixture model. We further compare the effect of using D andDMS as the anatomical connectivity estimates. Fig. 3.6 shows the numberof detected ROIs for different p-value thresholds. Our approach is denotedas RW for the first iteration where the full anatomical adjacency matrix isused, and RWit for the following iterations with partial adjacency matrices.For the same specificity, our approach provided higher detection sensitiv-ity than using the activation likelihoods given by the GGG mixture model,implicating the advantage of incorporating AC into activation detection.Iterating the procedure with refined AC estimates considerably improveddetection. The procedure was iterated 100 times and for clarity, only themean and standard deviations for the iterations with partial adjacency ma-trices were provided. It was observed that group activation map stabilizedafter 100 iterations, with only a couple of parcels changing labels in furtheriterations. This additional improvement provided by the refinement of theAC estimate suggests that only a subset of fibers is employed during a giventask. Hence, isolating the utilized fibers can be beneficial. The improvementcould also be partly due to removal of false fiber tracts arising from trac-tography errors. GLM (FDR corrected) provides more detection than ourmethod at more liberal thresholds, but its performance varies considerablywith p-value thresholds. In contrast, our method provides more consistentresults across p-values. Comparing the results of using the original anatomi-cal adjacency matrix versus the multi-step counterpart (denoted as RW-MSfor the first, RW-MSit for the following iterations) suggests that some of theartifactual connections generated by the multi-step approach likely do notpertain to task activation.Qualitatively, the ROIs detected across all experimental conditions inhealthy controls largely match regions known to be involved in context mem-ory tasks [139], which further validates our method. As shown in Fig. 3.7,the areas detected across all conditions include superior frontal gyrus, middlefrontal gyrus, inferior frontal gyrus, supramarginal gyrus, precentral gyrus,603.2. Modelling Anatomical Connectivity0 0.01 0.02 0.03 0.04 0.05020406080100120140  GGGGLMRWµ(RWit)σ(RWit)RW−MSµ(RW−MSit)σ(RW−MSit)Number of active parcelsp−value threshold (corrected)Figure 3.6: Quantitative results of the proposed multimodal activation de-tection technique.postcentral gyrus, angular gyrus, lateral occipital cortex, occipital pole, cin-gulate gyrus, lingual gyrus, hippocampus and insula.To assess the significant activation differences between schizophrenia pa-tients and control subjects, we compared two different types of context mem-ory recall: (1) self-other source information monitoring (did I say this wordor did I hear it?) and (2) task information monitoring (did I produce asemantic associate of this word or did I read it?). Both types of contextmemory are thought to be impaired in schizophrenia. This comparison cor-responds to the contrast between associating/reading and solving/hearingconditions. We employed a max-t permutation test [67] on the pA values ofthe two groups for this contrast and observed a significant difference (p<0.05,corrected) in the left hippocampus. Activity in this region was higher in thesource monitoring condition for schizophrenia patients relative to controlsbut lower in the task monitoring condition relative to controls. Applyingthe same test on the results of GLM or GGG mixture modelling approachfailed to detect any group differences up to p<0.2. The left hippocampusis involved in reactivation and association of stored semantic knowledge toconsolidate new information into existing semantic frameworks [140–142].The implication of the significant difference in the left hippocampus is that,aberrant activity in this region could lead to different manifestations of poor613.2. Modelling Anatomical ConnectivityFigure 3.7: Qualitative results of the proposed multimodal activation detec-tion method at p-value < 0.05. Green=detected by RW only. Blue=detectedby GGG only. Red=detected both by RW and GGG. Each row correspondsto an experimental condition, top-to-bottom: associating, hearing, solving,reading.623.3. Summaryperformance in context memory. Reduced activity during task monitoringcould relate to the episodic memory impairments commonly observed inschizophrenia [143], which implies that underactivity in this region wouldlead to reduced reactivation of stored semantic knowledge for context mem-ory. Increased activity during source monitoring could relate to the sourcememory impairments commonly observed in schizophrenia [144], such thatoveractivity in this region would lead to increased perceptual context forboth self and other source information, leading to more difficulty in distin-guishing between these two sources in context memory.3.2.7 DiscussionIn this section, we proposed a novel fiber connectivity integrated approachfor group activation inference. On real data of healthy subjects, we demon-strated that integrating dMRI and fMRI significantly increases sensitivityin detecting group activation compared to analyzing fMRI data alone. Wefurther showed that incorporating a refined connectome comprising anatom-ical connections linked only to the estimated active brain regions results inimproved performance over using the full estimated connectome. Finally, wepresented novel findings in activation differences between healthy controlsand schizophrenia patients that were missed with standard methods. Ourmultimodal integration strategy thus holds great promise for brain activityanalysis.3.3 SummaryIn this chapter, we proposed two multimodal integration methods to improvethe sensitivity in activation detection from fMRI data. The first method in-volves the fusion of RS fMRI and task fMRI data, and aims at deconfoundingthe effects of RS modulations from task fMRI data. Our results imply thatRS activity also contributes to the noise seen in task fMRI data, and thatits effects can be diminished by considering RS fMRI and task fMRI jointlyin a multimodal framework. Our second proposed method involves the fu-sion of dMRI and task fMRI data. We incorporated AC information intotask activation detection motivated by the fact that fiber pathways serveas the physical substrate for functional interactions, and showed that thismultimodal approach increases sensitivity in activation detection. Given thestrong noise in fMRI data, we illustrated that multimodal methods providepromising solutions to regularize the challenging problem of brain activationdetection.63Chapter 4Multimodal BrainSubnetwork IdentificationThe brain is known to consist of spatially disparate regions that are func-tionally connected as subnetworks [145]. Unraveling this modular structureof the brain is a challenging endeavor. To identify subnetworks from fMRIdata, seed-based correlation [22] and ICA [86] are primarily used. Morerecently, graph theoretical approaches [146] are gaining popularity as theyfacilitate compact brain representation and enable the reduction of complexbrain interactions into intuitive summary network measures. Under thisformalism, the brain is modelled as a graph, where brain regions and theirconnections correspond to nodes and edges, respectively.As described in Section 1.6.2, AC and FC estimates suffer from opposingartifacts. Namely, FC estimates tend to contain many false positive connec-tions due to several confounding factors involved, whereas AC estimates of-ten suffer from false negatives due to premature termination of reconstructedtracts at fiber crossings. We thus hypothesize that multimodal integrationfor subnetwork identification should prove beneficial since incorporating ACwould reduce the effects of noise-induced FC and vice versa.In this chapter, we propose a technique based on replicator dynamics(RD) that facilitates integration of fMRI and dMRI data in identifying brainsubnetworks. We start by reviewing the properties of RD [147], which onlyChapter 4 is adapted from the following papers:B. Yoldemir, B. Ng, R. Abugharbieh. Overlapping replicator dynamics for functionalsubnetwork identification. In Proceedings of the International Conference on MedicalImage Computing and Computer Assisted Intervention (MICCAI’13), LNCS, vol. 8150,pp. 682-689, Nagoya, Japan, September 2013.B. Yoldemir, B. Ng, R. Abugharbieh. Coupled stable overlapping replicator dynamicsfor multimodal brain subnetwork identification. In Proceedings of the 24th Biennial In-ternational Conference on Information Processing in Medical Imaging (IPMI’15), LNCS,vol. 9123, pp. 770-781, Isle of Skye, Scotland, June 2015.B. Yoldemir, B. Ng, R. Abugharbieh. Stable overlapping replicator dynamics for braincommunity detection. IEEE Transactions on Medical Imaging, vol. 35, no. 2, pp. 529-538,February 2016.644.1. Replicator Dynamicspermits unimodal subnetwork identification (Section 4.1). We then describea sex-differentiated formulation of RD [148] (Section 4.2) that permits FCand AC to be jointly modelled. In particular, this formulation extracts cou-pled subnetwork pairs based on fitness measures of both females and males,which in the present context, correspond to FC and AC estimates. We referto this formulation as coupled RD (CRD). Importantly, we prove that fordense graphs, CRD always selects the same set of nodes for each pair ofidentified subnetworks (Section 4.2). Hence, the effective output of CRDis actually a single set of subnetworks with FC and AC integrated. CRDalone is prone to noise and tends to falsely divide the subnetworks. Wethus adopt a graph incrementation scheme to facilitate merging of subnet-work components that might be falsely split and present a procedure forincorporating stability selection [149] to statistically control for inclusion offalse nodes arising from noise and overmerging (Section 4.3). We refer tothis extension of CRD as coupled stable RD (CSRD). Further, we describehow a graph augmentation strategy can be integrated into CSRD to enableoverlapping subnetwork identification (Section 4.4). In aggregate, the re-sulting technique: (i) facilitates multimodal integration, (ii) can operate onweighted graphs, (iii) allows for subnetwork overlaps, (iv) has an intrinsiccriterion for determining the number of subnetworks, (v) does not force allnodes to be part of a subnetwork, and (vi) statistically controls for falsenode inclusion. We refer to this technique as coupled stable overlapping RD(CSORD).4.1 Replicator DynamicsRD is a concept that originated from theoretical biology for modelling theevolution of interacting and self-replicating entities [147]. Let w(k) =(w1(k), · · · ,wd(k))T be a d × 1 vector, where wj(k) is the proportion ofallele that are of type j in the gene pool during generation k, and d is thenumber of alleles. Further, let C be a d× d matrix with each element, Cij ,reflecting the fitness of a genotype, i.e. a pair of alleles i and j. Under theassumption of natural selection, w(k+1) is given by the replicator equation[147]:w(k + 1) =w(k) ◦Cw(k)w(k)TCw(k), (4.1)where ◦ denotes element-wise multiplication. Eq. (4.1) is shown to solvethe following non-convex optimization problem [129]:maxwwTCw s.t. ||w||1 = 1,w ≥ 0. (4.2)654.2. Coupled Replicator DynamicsThe fundamental theorem of natural selection states that if C is real-valued, symmetric, and non-negative, w(k)TCw(k) strictly increases with kalong any nonstationary trajectory w(k) until it converges to a strict localmaximum as (4.1) is iterated [150]. In the more general case of any real-valued C, local maxima of w(k)TCw(k) are guaranteed to attract nearbytrajectories w(k) asymptotically as (4.1) is iterated [151]. wj(k + 1) willequal wj(k) upon convergence thus (Cw(k))j = w(k)TCw(k) for all j,where (·)j denotes the jth element of the corresponding vector. In graph-theoretic terms with C being a similarity matrix, this implies that all se-lected nodes (i.e. nodes with non-zero w(k)) will have the same weightedaverage correlation with one another. Hence, the grouping of the nodes isbased on mutual similarity among the nodes. With w constrained to lieon the standard simplex, i.e. ||w||1 = 1, w ≥ 0, sparse w is encouragedsince restricting the l1 norm induces sparsity [112]. wj(k) > 0 hence indi-cates that allele j persisted in generation k. The local maximum to whichRD converges highly depends on w(0). Since the vicinity at which the lo-cal maxima reside is unknown a priori, an unbiased way of setting w(0)is to assign all of its elements to 1/d, which tends to find the communitywith the highest mutual similarity among nodes [129]. Other subnetworkscan be found by reapplying (4.1) after removing the nodes in the identifiedsubnetworks [152].4.2 Coupled Replicator DynamicsSeparating C into d× d genotype fitness matrices F for females and M formales to model sex differences, the survival probabilities of the alleles inthe female gene pool, p(k), and in male gene pool, q(k), based on naturalselection can be estimated using the following multiplicative updates [148]:p(k + 1) = 0.5p(k) ◦ Fq(k) + q(k) ◦ Fp(k)p(k)TFq(k)(4.3)q(k + 1) = 0.5q(k) ◦Mp(k) + p(k) ◦Mq(k)p(k)TMq(k). (4.4)Analogous to RD, F and M are assumed to be real, symmetric, andnon-negative. The non-negativity constraint enforces p(k) and q(k) to be664.2. Coupled Replicator Dynamicsgreater than or equal to 0. Also,∑j pj(k) = 1 since:0.5 ·∑j(pj(k)(Fq(k))j + qj(k)(Fp(k))j) =0.5 · (p(k)TFq(k) + q(k)TFp(k)) =0.5 · (p(k)TFq(k) + p(k)TFq(k)) =p(k)TFq(k),(4.5)where we have exerted the property that a scalar and its transpose areequal, i.e. q(k)TFp(k) = p(k)TFq(k). Similarly,∑j qj(k) = 1. Thus, bymatching the denominator of (4.3) and (4.4) with that of (4.1), one cansee that (4.3) and (4.4) together form a local maximizer of the followingoptimization problem:maxp,qpTFq + pTMq s.t. ||p||1 = 1,p ≥ 0, ||q||1 = 1,q ≥ 0, (4.6)where the l1 constraint on p and q enforces sparsity. Since p is multiplied toboth F and M in (4.6), the optimal p would jointly account for the genotypefitness of both sexes, and the same goes for q. This coupling property ofCRD can be more clearly seen by examining (4.3) and (4.4). Specifically,information in M is propagated to p through q, and vice versa. Further,we highlight that the sparsity patterns of p and q tend to perfectly matchupon convergence for dense graphs. Specifically, assume pi(k) = 0 uponconvergence. Based on (4.3), this implies pi(k)◦Fqi(k)+qi(k)◦Fpi(k) = 0.pi(k) ◦ Fqi(k) must equal 0 since pi(k) = 0, thus qi(k) ◦ Fpi(k) mustalso equal 0. This condition is satisfied only if qi(k) = 0 or Fpi(k) = 0.Fpi(k) = 0 requires Fij to be exactly zero for each surviving allele j. Inboth our synthetic and real data, F is never sparse enough to warrant thiscondition. Hence, qi(k) ◦ Fpi(k) = 0 can only arise by having qi(k) = 0,which enforces p and q to have the same sparsity patterns. The nodes in theidentified subnetwork are thus given by the nonzero elements of either p orq. Again, other subnetworks can be found by reapplying (4.3) and (4.4) afterremoving the nodes of the identified subnetworks from the graphs. However,this strategy modifies the original graph structure, thereby affecting theremaining local maxima. We describe a more principled way of identifyingthe other subnetworks that permits subnetwork overlaps and preserves thelocation of the remaining local maxima in Section 4.4.674.3. Coupled Stable Replicator Dynamics4.3 Coupled Stable Replicator Dynamics4.3.1 Graph IncrementationSparsity of the subnetworks identified by CRD is intrinsically determinedbased on the mutual similarity among nodes. This property is desirable un-der the ideal case of noiseless graphs. In practice, such a strict criterion forgrouping nodes could easily result in a subnetwork being falsely split intocomponents [153], since even small perturbations to F and M would rendercertain nodes non-mutually connected. To reduce this effect, one strategyis to add a constant, η, to the off-diagonal elements of F and M. Usingan extreme example for intuition, if we add 10,000 to F and M with valuesranging from 0 to 1, the original differences between the elements of F (andM) would be negligible. Applying CRD to such matrices would result inall nodes being assigned to a single subnetwork, i.e. all elements of p andq would attain nonzero values. If we instead use a η of magnitude similarto the values in F and M, only the smaller differences in F and M wouldbe negligible. This strategy thus artificially increases the mutual similar-ity among all nodes and allows adjusting the sparsity level of the identifiedsubnetworks by changing η. Choosing the optimal η is non-trivial. Also,noise in F and M could result in false nodes being included in the subnet-works. To jointly deal with these problems, we present next a procedure forintegrating stability selection into CRD.4.3.2 Stability SelectionThe idea behind stability selection is that if we bootstrap the data manytimes and perform subnetwork identification on each bootstrap sample,nodes that truly belong to the same subnetwork are likely to be jointlyselected over a large fraction of bootstrap samples, whereas false nodes areunlikely to be persistently selected [149]. To incorporate stability selectionin refining each subnetwork identified by CRD, we first generate a set ofmatrices {Fη} by adding η to the off-diagonal elements of F. A range of ηfrom 0 to d ·maxi,j |Fij | at a step size of κ is used, where κ is chosen smallenough to ensure that no more than one node is added to a subnetwork ateach increment of η. The same procedure is applied to M to generate {Mη}.We then apply CRD to each (Fη,Mη) pair with p(0) and q(0) set to theCRD solution of the previous η increment (with 1/d added to all elementsand renormalized to sum to 1), which enforces CRD to converge to a localmaximum in the vicinity of the previous η increment for retaining subnet-work correspondence. We note that adding 1/d is critical since zeros in p(k)684.4. Coupled Overlapping Stable Replicator Dynamicsand q(k) remain as zeros as (4.3) and (4.4) are iterated, which prohibits newnodes from being added. We discard subnetworks that contain more than10% of the nodes assuming a subnetwork would not span more than 10%of the brain. Next, we generate 100 bootstrap samples of F and M (as ex-plained in Section 4.5), and compute{Fbη}and{Mbη}from each bootstrapsample Fb and Mb. We then apply CRD on each (Fbη,Mbη) pair. For each η,we estimate the selection probability of each node as the proportion of boot-strap samples in which the given node attains a nonzero weight [149]. Thisestimation can be done based on either p or q since their sparsity patternsperfectly match. Finally, we threshold the selection probabilities to identifynodes that belong to the same subnetwork. A threshold τ that bounds theexpected number of false nodes in a subnetwork, E, is given by [149]τ ≤(q2Ed+ 1)/2, (4.7)where q is the average number of nodes per subnetwork, which can be es-timated as the sum of selection probabilities of all nodes averaged over η.We set E to 1 so that the average number of false nodes included in a sub-network is statistically controlled to be less than or equal to 1. We declarethe set of nodes with selection probabilities higher than the resulting τ forany η to be a subnetwork [149]. We note that each subnetwork identified byCRD is refined independently using the above procedure.4.4 Coupled Overlapping Stable ReplicatorDynamicsSince CRD is deterministic, reapplying it with the same initialization willconverge to the same local maximum. Here, we adopt a graph augmenta-tion strategy that enables destabilization of previously-found local maximawithout altering the location of the others. Reapplying CRD on this aug-mented graph would thus converge to one of the remaining local maxima.694.5. Parameter Selection and ImplementationThe required graph augmentation is as follows [151]:Faugij =Fij if i, j ≤ dα if j > d and i /∈ Sj−dβ if i, j > d and i = jγij if i > d and j ∈ Si−d0 otherwiseγij = σij1|Si−d|∑m∈Si−dFmj , σij > 1, (4.8)where α > β, β = maxi 6=j Fij , Sl is the set of nodes of the lth identified sub-network, and | · | denotes cardinality. M is similarly augmented to generateMaug. In effect, this is equivalent to adding artificial nodes to the graphsby appending F and M with new rows and new columns, and extendingweighted edges to the original nodes such that the previously-found localmaxima provably become unstable. Given S1 found using CRD and refinedby stability selection, finding S2 proceeds by adding a new row and a newcolumn to F and M based on (4.8), applying CRD on the new (Faug,Maug)pair, and refining the resulting subnetwork using stability selection. Remain-ing Sl can be similarly extracted. We highlight that an intrinsic criterionfor terminating further subnetwork extraction is to stop if after convergence,p(k)TFq(k)+p(k)TMq(k) ≤ p(0)TFq(0)+p(0)TMq(0), since this suggeststhat no further solutions of (4.6) are present. In practice, we found on syn-thetic data that using a·(p(0)TFq(0)+p(0)TMq(0)), a > 1, as the stoppingcriterion is more robust to noisy F and M. A rigorous way of choosing a isdiscussed in Section Parameter Selection and ImplementationIn the context of brain subnetwork identification, F and M correspond toestimates of FC and AC, respectively, and alleles correspond to brain ROIs.Our goal is to find subnetworks of ROIs that are highly inter-connected bothfunctionally and anatomically. For F, we use the conventional Pearson’scorrelation between average ROI time series. To ensure F is non-negative asrequired for the properties of CRD to hold, we set negative Fij to zero giventhe currently unclear interpretation of negative FC. For M, we use the fibercount between ROIs following our discussion in Section 2.3.2. The maindiagonals of F and M are set to zero to avoid self-connections [129]. Given704.6. Materialsthe substantially larger magnitudes in M compared to F, the subnetworkdetection process would be dominated by M. To mitigate this effect, welinearly rescale M such that the values in F and M are in the same range.An open problem in subnetwork identification is the choice on the num-ber of subnetworks. Although p(k)TFq(k) + p(k)TMq(k) ≤ p(0)TFq(0) +p(0)TMq(0) is theoretically justified, noisy F and M could lead to spuri-ous subnetworks with p(k)TFq(k) + p(k)TMq(k) being slightly above thethreshold. Hence, there is a need for a stricter threshold, a · (p(0)TFq(0) +p(0)TMq(0)), a ≥ 1. To set a, we generate 500 synthetic datasets andextract subnetworks over a range of a from 1 to 10. The value of a thatminimizes the difference between the number of estimated subnetworks andthe ground truth is declared as the optimal a. We note that the syntheticdatasets used for selecting a are separate from those used for evaluating theperformance of CSORD.To generate bootstrap samples of F for applying stability selection, weadopt a parametric bootstrap approach [154] since the spatiotemporal corre-lation structure of fMRI data would not be retained with regular bootstrap,where time points are independently sampled with replacement. Let X bea t × d fMRI time series matrix, where t is the number of time samples.We randomly draw 100 Xb of size t × d from a multivariate normal distri-bution with zero mean and covariance, Fpd, where Fpd is a positive definiteapproximation of F, computed by setting non-positive eigenvalues of F to10−16. Fb is then computed as the Pearson’s correlation of Xb. As for M,one approach for generating bootstrap samples is to randomly select diffu-sion gradient directions with replacement to generate bootstrapped dMRIvolumes and reapply tractography [155], but this is very computationallyintensive. Instead, we use the same parametric bootstrap procedure to gen-erate bootstrap samples of M.4.6 Materials4.6.1 Synthetic DataWe generated 500 synthetic datasets that cover a wide variety of networkconfigurations. Each dataset comprised d = 200 regions and 4 scans of 1200time points as in the real data. We set the number of subnetworks, N , ineach dataset to a random value between 10 and 20. The number of ROIsin each subnetwork was set to dd/Ne + c with c being a random numberbetween 1 and 5, which ensures that the total number of ROIs across sub-networks is greater than d, hence guaranteeing the presence of subnetwork714.7. Resultsoverlaps. ROIs were randomly assigned to subnetworks, and the same ROIwas allowed to repeat across subnetworks. We opted to used random net-work configurations since the ground truth number of subnetworks, numberof nodes within each subnetwork, and amount of subnetwork overlap areunknown. Our strategy is thus to try to exhaust the possibilities. For eachnetwork configuration, we generated a d×d adjacency matrix, Σ, which wastaken as the ground truth. Next, we built a noise-free AC matrix, Σa, byrandomly setting p1% of the values in Σ to 0 to model how AC estimatesare prone to false negatives, and built a noise-free FC matrix, Σf , by ran-domly setting p2% of the values in Σ to 1 to model how FC estimates areprone to false positives. p1 and p2 were randomly chosen from [0, 20]. Twosets of time series were then generated by drawing random samples fromN (0,Σa) and N (0,Σf ). Finally, AC and FC matrices were simulated byadding Gaussian noise to the time series with SNR randomly set between-6 and -3 dB, and computing the Pearson’s correlation of these noisy timeseries. We chose the typical SNR levels seen in task fMRI data [23, 24] sinceSNR of RS fMRI data is hard to determine.4.6.2 Real DatadMRI, RS fMRI and task fMRI data from the HCP dataset were used in theexperiments [123]. Details and preprocessing of the data and the generationof the brain parcellation map are presented in Sections 2.2.3 and Results4.7.1 Synthetic Data ResultsWe compared CSORD against CSRD, CRD, stable RD (SRD), RD, NC[156], and multi-view spectral clustering (MVSC) [52]. SRD is analogous toCSRD, except (4.1) is used instead of (4.3) and (4.4) in finding subnetworks.Among the contrasted techniques, CSORD, CSRD, CRD and MVSC per-mit multimodal subnetwork identification, whereas SRD, RD, and NC arelimited to a single modality. For these unimodal techniques, we focused ontheir performance in extracting subnetworks from FC estimates. On 500synthetic datasets, we assessed each contrasted technique by computing theDice coefficient (DC) between the estimated and ground truth subnetworksasDC = 2|Sest ∩ Sgnd||Sest|+ |Sgnd| , (4.9)724.7. Resultswhere Sest is the set of ROIs in an estimated subnetwork, and Sgnd is the setof ROIs in the corresponding ground truth subnetwork matched to Sest usingHungarian clustering [157]. DC of unmatched subnetworks was set to zero.CSORD achieved significantly higher DC than each contrasted technique atp< 10−40 based on Wilcoxon signed rank test (Fig. 4.1).CSORD CSRD CRD SRD RD NC MVSC00. DCFigure 4.1: Subnetwork identification accuracy on synthetic data.4.7.2 Real Data ResultsSince ground truth subnetworks are unknown for real data, we based ourvalidation on test-retest reliability and task classification accuracy. We es-timated test-retest reliability as the intra-class correlation (ICC) over fourgraph metrics. Specifically, we first computed the clustering coefficient,transitivity, global efficiency, and local efficiency [158] for each subnetworkextracted from RS fMRI data (and dMRI data for multimodal methods) ofeach scan session of a given subject. We then averaged the values of a givenmetric across each subject’s subnetworks for each of the two scan sessions.Collecting these average values of each metric across subjects, we computedthe ICC between the two scan sessions. We finally estimated test-retestreliability as the average of resulting four ICC values.High test-retest reliability is a desired property, but could arise fromcommon noise between scan sessions. We thus additionally assessed thetask classification accuracy of the contrasted techniques as follows. Wefirst extracted group subnetworks from RS fMRI and dMRI data. We thencomputed F of each task from each subject’s task fMRI dataset, and scaled734.7. Resultsits elements asFsij = Fij ·|mi ∩mj ||mi ∪mj | , (4.10)where mi and mj are the group subnetwork memberships of nodes i andj based on RS fMRI and dMRI data. This scaling weighs down Fij ifnodes i and j were assigned to different group subnetworks, which in effect,propagates the group subnetwork structure onto F. The weighted degreeof each brain region computed from Fs was taken as a feature, resulting ina d × 1 feature vector for each task. A support vector machine was usedfor classification. For estimating classification accuracy, we subsampled thedata over 100 random splits with 35 subjects used for training and 5 fortesting.We adopted the above scaling scheme of F, which was originally usedfor extending the definition of modularity to overlapping subnetworks [159],since the notion of within and between subnetwork connections become fuzzywhen subnetworks are allowed to overlap. Also, this evaluation schemeensures the data used for subnetwork identification (RS fMRI and dMRIdata) are independent from the data used for classification performanceassessment (task fMRI data). The assumption is that tasks induce onlysmall functional changes to intrinsic brain connectivity. This hypothesisis supported by the strong resemblance between RS and task subnetworks[78], as well as how energy consumption in the brain is mainly for supportingongoing activity with task-evoked response constituting less than 5% of thetotal [76].Results on test-retest reliability and task classification accuracy are shownin Fig. 4.2. Comparing RD’s and CRD’s performance suggests that simplyincorporating AC without considering the limitations of RD produces onlymarginal benefits. In fact, a much larger gain was obtained by controllingfor false node inclusions with SRD. Nevertheless, by jointly incorporatingAC and controlling for false node inclusion as well as enabling subnetworkoverlaps, CSORD achieved the highest task classification accuracy and ison par with the best contrasted techniques in terms of test-retest reliability.We note that average FC has been argued to contain no dynamical informa-tion. It is true that average connectivity matrix does not capture dynamicalconnectivity, but the overlapping structure of subnetworks is well reflectedby the relative magnitude of its elements, and can be readily extracted asshown by our results.Qualitatively, CSORD identified all commonly found subnetworks in theliterature [160]. The extracted visuospatial subnetwork and right execu-tive control subnetwork are shown in Figs. 4.3 and 4.4 as exemplar re-744.7. ResultsAverage9ICCTask9classification9accuracy00. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5MVSCNCRD CRDCSRDSRDCSORDFigure 4.2: Subnetwork identification results on real data. CSORD achievedthe best classification accuracy and high test-retest reliability.sults. CSORD also identified a subnetwork (Fig. 4.5) that comprises thebasal ganglia, thalamus, higher visual cortex, motor areas, and part of thecerebellum. These regions constitute the visual corticostriatal loop [161],striato-thalamo-cortical (STC) loop, and cerebello-thalamo-cortical (CTC)loop [162], which are rarely identified as a single subnetwork. Instead, state-of-the-art techniques typically declare the basal ganglia and cerebellum asseparate subnetworks [78, 163], failing to capture cortical projections fromsubcortical regions.Figure 4.3: Visuospatial subnetwork identified by CSORD at the group level.754.7. ResultsFigure 4.4: Right executive control subnetwork identified by CSORD at thegroup level.Figure 4.5: A group subnetwork identified by CSORD that comprises thevisual corticostriatal, CTC, and STC loops. Note that the cerebellum isoutside the presented views.764.8. Discussion4.8 DiscussionWe proposed CSORD to combine information from fMRI and dMRI datain identifying overlapping brain subnetworks. On synthetic data coveringa diverse set of network configurations, we showed that CSORD providessignificantly higher subnetwork identification accuracy than a number ofstate-of-the-art techniques. On real data, we illustrated that CSORD is notonly more robust to inter-session variability than most of the contrastedtechniques, but it also achieved the highest classification accuracy. Impor-tantly, standard techniques typically declare the basal ganglia, cerebellum,and motor cortex as separate subnetworks, despite their known joint in-volvements in various motor loops. In contrast, CSORD was able to extracta single subnetwork comprising these well-known interacting motor regions.Collectively, our results show that multimodal integration is a promisingapproach for brain subnetwork identification.77Chapter 5ConclusionsWe presented in this thesis our technical contributions for the assessmentof AC, and joint analysis of AC and FC for investigating functional segre-gation and integration in the human brain. Our techniques were based onmultimodal integration of RS fMRI, task fMRI and dMRI data. We con-clude in this section with a discussion of the strengths and weaknesses ofour methods and some promising directions to explore.In Chapter 2, we proposed the use of a super-resolution technique indMRI to ameliorate the adverse effects of low spatial resolution on AC es-timation. We showed that our approach qualitatively improves the fibertracts and tract-density maps generated and significantly improves the cor-relation between AC and FC estimates on real data. We further showedthat the traditional way of estimating AC from dMRI data, i.e. streamlinetractography on diffusion tensors, leads to understating the AC-FC consis-tency, and using ODFs over diffusion tensors and global tractography overstreamline tractography help increase the consistency between AC and FCestimates. Finally, we illustrated the effect of the choice of AC metric on theAC-FC consistency and showed that fiber count provides a good compromisebetween high AC-FC consistency and low fiber length bias.In Chapter 3, we addressed the challenging problem of brain activationdetection and proposed two multimodal strategies to increase sensitivity indetecting task activation. Our first method involved the fusion of RS fMRIand task fMRI data with the aim of deconfounding the effects of RS activityfrom task fMRI data. We showed that the removal of estimated RS modu-lations from task fMRI data significantly improves activation detection. Weused a linear regression model for the removal of spontaneous fluctuationsmotivated by the work of Fox et al. [75], where the results imply that task-evoked signal and RS activity are roughly linearly superimposed in the brain.Our second proposed method involved the fusion of dMRI and task fMRIdata. We showed that our fiber connectivity integrated activation detectionapproach increases sensitivity in detecting activation in controls, and alsoenables the detection of significant group differences between controls andschizophrenia patients. To the best of our knowledge, this is the first work785.1. Future Workin the literature to regularize activation detection using AC information.In Chapter 4, we proposed CSORD to identify groups of brain regionsthat exhibit high inter-connectivity both functionally and anatomically. Ourproposed method draws on advances from theoretical biology [147], evolu-tionary game theory [151] and statistics [149]. There are several advantagesof the proposed method: (i) it facilitates multimodal integration, (ii) it canoperate on weighted graphs, (iii) it allows for subnetwork overlaps, (iv) ithas an intrinsic criterion for determining the number of subnetworks, (v) itdoes not force all nodes to be part of a subnetwork, and (vi) it statisticallycontrols for false node inclusion. We demonstrated how this method enablesmore accurate detection of subnetworks on synthetic data compared to anumber of state-of-the-art approaches. Also, we showed that CSORD iden-tifies more reproducible subnetworks and facilitates successful classificationof the tasks subjects are involved in by analyzing the subnetworks recruited.5.1 Future WorkThroughout Chapter 2, we used Pearson’s correlation to assess the consis-tency between AC and FC. This linear model was used instead of higher-order nonlinear models in estimating AC-FC consistency for easier inter-pretability, and based on the empirical finding that AC and FC in thehuman brain are roughly linearly related [64, 97] when AC and FC arequantified using the metrics we have adopted in this thesis. Nonetheless, un-derstanding the exact nature of this relationship warrants further research.It should be emphasized that there is currently no well-accepted definitionof “anatomical connectivity strength”. All AC metrics we have consideredin this thesis constitute indirect measures based on tractography that cap-ture different aspects of the connectivity pattern. Morphological measuressuch as cortical thickness computed from structural MRI data have alsobeen employed for this purpose [164], however these are also proxies for ACstrength. Also, cortical thickness cannot be computed for subcortical ar-eas unlike tractography-based metrics. Devising more direct AC metrics iscurrently an open research problem beyond the scope of this thesis.In Chapter 3, we used a linear model to deconfound RS modulationsfrom task fMRI data. Our work was motivated by prior work which impliesa largely linear superposition of task-evoked signal and RS modulations inthe brain [75]. Our results in Section 3.1.4 support this linearity claim. How-ever, we believe that understanding whether there are additional nonlinearinteractions between the task-evoked signal and RS activity warrants further795.1. Future Workresearch. Sensitivity in brain activation detection can then presumably befurther increased by modelling such additional interactions.A limitation of CSORD presented in Chapter 4 is its high computationalcomplexity. The method involves running an iterative scheme for each boot-strap sample, thus may not easily extend to voxel-level analysis. Neverthe-less, the refinement of each subnetwork using stability selection is indepen-dently performed, hence the method may benefit from parallel programmingto speed up the computation. Identifying subject-specific subnetworks usingCSORD and examining whether the topology of these subnetworks can belinked to e.g. disease severity or genetic factors is a promising direction toexplore in the future.80Bibliography[1] J. Hooper and D. Teresi. The three-pound universe. Macmillan, 1986.[2] Mapping connections: An understanding of neurological conditions inCanada. Public Health Agency of Canada, 2014.[3] The global burden of disease: Generating evidence, guiding policy.Institute for Health Metrics and Evaluation, University of Washington,2012.[4] V. Hesselmann, B. Sorger, R. Girnus, K. Lasek, M. Maarouf,C. Wedekind, J. Bunke, O. Schulte, B. Krug, K. Lackner, andV. Sturm. Intraoperative functional MRI as a new approach to moni-tor deep brain stimulation in Parkinson’s disease. European Radiology,14(4):686–690, 2004.[5] G.D. Honey, E. Pomarol-Clotet, P.R. Corlett, R.A.E. Honey, P.J.Mckenna, E.T. Bullmore, and P.C. Fletcher. Functional dysconnectiv-ity in schizophrenia associated with attentional modulation of motorfunction. Brain, 128(11):2597–2611, 2005.[6] T.C. Baghai, H.-J. Moller, and R. Rupprecht. Recent progress inpharmacological and non-pharmacological treatment options of majordepression. Current Pharmaceutical Design, 12(4):503–515, 2006.[7] S.Y. Bookheimer, M.H. Strojwas, M.S. Cohen, A.M. Saunders, M.A.Pericak-Vance, J.C. Mazziotta, and G.W. Small. Patterns of brainactivation in people at risk for Alzheimer’s disease. New EnglandJournal of Medicine, 343(7):450–456, 2000.[8] M.D. Greicius, K. Supekar, V. Menon, and R.F. Dougherty. Resting-state functional connectivity reflects structural connectivity in the de-fault mode network. Cerebral Cortex, 19(1):72–78, 2009.[9] S. Jbabdi and H. Johansen-Berg. Tractography: Where do we go fromhere? Brain Connectivity, 1(3):169–183, 2011.81Bibliography[10] K.J. Friston. Models of brain function in neuroimaging. Annual Reviewof Psychology, 56(1):57–87, 2005.[11] G. de Marco. Effective connectivity and brain modeling by fMRI.Advanced Studies in Biology, 1(3):139–144, 2009.[12] O. Sporns, G. Tononi, and R. Ko¨tter. The human connectome: astructural description of the human brain. PLoS Computational Biol-ogy, 1(4):e42, 2005.[13] R.C. Craddock, S. Jbabdi, C.-G. Yan, J.T. Vogelstein, F.X. Castel-lanos, A. Di Martino, C. Kelly, K. Heberlein, S. Colcombe, and M.P.Milham. Imaging human connectomes at the macroscale. NatureMethods, 10(6):524–539, 2013.[14] P. J. Basser, J. Mattiello, and D. Le Bihan. Estimation of the effectiveself-diffusion tensor from the NMR spin echo. Journal of MagneticResonance, Series B, 103(3):247–254.[15] H.-E. Assemlal, D. Tschumperle´, L. Brun, and K. Siddiqi. Recent ad-vances in diffusion MRI modeling: Angular and radial reconstruction.Medical Image Analysis, 15(4):369–396, 2011.[16] P. Hagmann, L. Jonasson, P. Maeder, J.-P. Thiran, V.J. Wedeen, andR. Meuli. Understanding diffusion MR imaging techniques: Fromscalar diffusion-weighted imaging to diffusion tensor imaging and be-yond. Radiographics, 26:S205–S223, 2006.[17] D. Le Bihan, C. Poupon, A. Amadon, and F. Lethimonnier. Artifactsand pitfalls in diffusion MRI. Journal of Magnetic Resonance Imaging,24(3):478–488, 2006.[18] J. Soares, P. Marques, V. Alves, and N. Sousa. A hitchhiker’s guideto diffusion tensor imaging. Frontiers in Neuroscience, 7(31), 2013.[19] S. Ogawa, T. M. Lee, A. R. Kay, and D. W. Tank. Brain magnetic res-onance imaging with contrast dependent on blood oxygenation. Pro-ceedings of the National Academy of Sciences, 87(24):9868–9872, 1990.[20] E. Amaro and G.J. Barker. Study design in fMRI: Basic principles.Brain and Cognition, 60(3):220–232, 2006.[21] D. Le Bihan, P. Jezzard, J. Haxby, N. Sadato, L. Rueckert, andV. Mattay. Functional magnetic resonance imaging of the brain. An-nals of Internal Medicine, 122(4):296–303, 1995.82Bibliography[22] B. Biswal, F. Zerrin Yetkin, V.M. Haughton, and J.S. Hyde. Func-tional connectivity in the motor cortex of resting human brain usingecho-planar MRI. Magnetic Resonance in Medicine, 34(4):537–541,1995.[23] G.E. Sarty. Computing brain activity maps from fMRI time-seriesimages. Cambridge University Press, 2007.[24] V. Calhoun. Preprocessing iii-iv: Anatomic transformation to stan-dard brain spaces. Accessed: 2015-08-18.[25] R.N.A. Henson, C. Buechel, O. Josephs, and K.J. Friston. The slice-timing problem in event-related fMRI. NeuroImage, 9:125, 1999.[26] K.J. Friston, J. Ashburner, C.D. Frith, J.-B. Poline, J.D. Heather, andR.S.J. Frackowiak. Spatial registration and normalization of images.Human Brain Mapping, 3(3):165–189, 1995.[27] A.P. Holmes, O. Josephs, C. Buchel, and K.J. Friston. Statisticalmodelling of low-frequency confounds in fMRI. NeuroImage, 5:S480,1997.[28] K.J. Worsley, C.H. Liao, J. Aston, V. Petre, G.H. Duncan, F. Morales,and A.C. Evans. A general statistical analysis for fMRI data. Neu-roImage, 15(1):1–15, 2002.[29] J.C. Mazziotta, A.W. Toga, A. Evans, P. Fox, and J. Lancaster. Aprobabilistic atlas of the human brain: Theory and rationale for its de-velopment: The international consortium for brain mapping (ICBM).NeuroImage, 2(2):89–101, 1995.[30] W.D. Penny, K.J. Friston, J.T. Ashburner, S.J. Kiebel, and T.E.Nichols. Statistical parametric mapping: the analysis of functionalbrain images. Academic Press, 2011.[31] S.M. Smith, M. Jenkinson, M.W. Woolrich, C.F. Beckmann, T.E.J.Behrens, H. Johansen-Berg, P.R. Bannister, M. De Luca, I. Drobn-jak, and D.E. Flitney. Advances in functional and structural MR im-age analysis and implementation as FSL. NeuroImage, 23:S208–S219,2004.83Bibliography[32] Y. Chao-Gan and Z. Yu-Feng. DPARSF: a MATLAB toolbox forpipeline data analysis of resting-state fMRI. Frontiers in SystemsNeuroscience, 4, 2010.[33] P.J. Basser and C. Pierpaoli. Microstructural and physiological fea-tures of tissues elucidated by quantitative-diffusion-tensor MRI. Jour-nal of Magnetic Resonance, 111(3):209–219, 1996.[34] P.J. Basser, J. Mattiello, and D. LeBihan. MR diffusion tensor spec-troscopy and imaging. Biophysical Journal, 66(1):259, 1994.[35] P.J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi. Invivo fiber tractography using DT-MRI data. Magnetic Resonance inMedicine, 44(4):625–632, 2000.[36] T.E. Conturo, N.F. Lori, T.S. Cull, E. Akbudak, A.Z. Snyder, J.S.Shimony, R.C. McKinstry, H. Burton, and M.E. Raichle. Trackingneuronal fiber pathways in the living human brain. Proceedings of theNational Academy of Sciences, 96(18):10422–10427, 1999.[37] D.K. Jones, A. Simmons, S.C.R. Williams, and M.A. Horsfield. Non-invasive assessment of axonal fiber connectivity in the human brain viadiffusion tensor MRI. Magnetic Resonance in Medicine, 42(1):37–41,1999.[38] S. Mori and J.-D. Tournier. Introduction to diffusion tensor imagingand higher order models. Academic Press, 2013.[39] P.T. Callaghan. NMR imaging, NMR diffraction and applications ofpulsed gradient spin echoes in porous media. Magnetic ResonanceImaging, 14(78):701–709, 1996.[40] D.S. Tuch. Q-ball imaging. Magnetic Resonance in Medicine,52(6):1358–1372, 2004.[41] I. Aganj, C. Lenglet, G. Sapiro, E. Yacoub, K. Ugurbil, and N. Harel.Reconstruction of the orientation distribution function in single- andmultiple-shell q-ball imaging within constant solid angle. MagneticResonance in Medicine, 64(2):554–566, 2010.[42] M. Reisert, I. Mader, C. Anastasopoulos, M. Weigel, S. Schnell, andV. Kiselev. Global fiber reconstruction becomes practical. NeuroIm-age, 54(2):955–962, 2011.84Bibliography[43] T.E.J. Behrens, H.J. Berg, S. Jbabdi, M.F.S. Rushworth, and M.W.Woolrich. Probabilistic diffusion tractography with multiple fibre ori-entations: What can we gain? NeuroImage, 34(1):144–155, 2007.[44] G.J.M. Parker, H.A. Haroon, and C.A.M. Wheeler-Kingshott. Aframework for a streamline-based probabilistic index of connectivity(PICo) using a structural interpretation of MRI diffusion measure-ments. Journal of Magnetic Resonance Imaging, 18(2):242–254, 2003.[45] D.K. Jones and C. Pierpaoli. Confidence mapping in diffusion tensormagnetic resonance imaging tractography using a bootstrap approach.Magnetic Resonance in Medicine, 53(5):1143–1149, 2005.[46] M. Xu, X. Zhang, Y. Wang, L. Ren, Z. Wen, Y. Xu, G. Gong, N. Xu,and H. Yang. Probabilistic brain fiber tractography on GPUs. InIEEE 26th International Parallel and Distributed Processing Sympo-sium Workshops, pages 742–751, 2012.[47] B. Jeurissen, A. Leemans, J.-D. Tournier, K. Jones, and J. Sijbers.Investigating the prevalence of complex fiber configurations in whitematter tissue with diffusion magnetic resonance imaging. HumanBrain Mapping, 34(11):2747–2766, 2013.[48] Marco Catani. From hodology to function. Brain, 130(3):602–605,2007.[49] Q. Cao, N. Shu, L. An, P. Wang, L. Sun, M.-R. Xia, J.-H. Wang,G.-L. Gong, Y.-F. Zang, and Y.-F. Wang. Probabilistic diffusiontractography and graph theory analysis reveal abnormal white mat-ter structural connectivity networks in drug-naive boys with atten-tion deficit/hyperactivity disorder. The Journal of Neuroscience,33(26):10676–10687, 2013.[50] P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C.J. Honey, V.J.Wedeen, and O. Sporns. Mapping the structural core of human cere-bral cortex. PLoS Biology, 6(7):e159, 07 2008.[51] V. Tomassini, S. Jbabdi, J.C. Klein, T.E.J. Behrens, C. Pozzilli,P.M. Matthews, M.F.S. Rushworth, and H. Johansen-Berg. Diffusion-weighted imaging tractography-based parcellation of the human lateralpremotor cortex identifies dorsal and ventral subregions with anatom-ical and functional specializations. The Journal of Neuroscience,27(38):10259–10269, 2007.85Bibliography[52] H. Chen, K. Li, D. Zhu, X. Jiang, Y. Yuan, P. Lv, T. Zhang, L. Guo,D. Shen, and T. Liu. Inferring group-wise consistent multimodal brainnetworks via multi-view spectral clustering. IEEE Transactions onMedical Imaging, 32(9):1576–1586, 2013.[53] D.K. Jones, T.R. Kno¨sche, and R. Turner. White matter integrity,fiber count, and other fallacies: the do’s and don’ts of diffusion MRI.NeuroImage, 73:239–254, 2013.[54] K.J. Friston, C.D. Frith, P.F. Liddle, and R.S.J. Frackowiak. Func-tional connectivity: the principal-component analysis of large PETdata sets. Journal of Cerebral Blood Flow and Metabolism, 13:5–14,1993.[55] M.D. Fox and M. Greicius. Clinical applications of resting state func-tional connectivity. Frontiers in Systems Neuroscience, 4, 2010.[56] G. Marrelec, A. Krainik, H. Duffau, M. Plgrini-Issac, S. Lehricy,J. Doyon, and H. Benali. Partial correlation for functional brain inter-activity investigation in functional MRI. NeuroImage, 32(1):228–237,2006.[57] F.T. Sun, L.M. Miller, and M. D’Esposito. Measuring interregionalfunctional connectivity using coherence and partial coherence analysesof fMRI data. NeuroImage, 21(2):647–658, 2004.[58] S.M. Smith, K.L. Miller, G. Salimi-Khorshidi, M. Webster, C.F. Beck-mann, T.E. Nichols, J.D. Ramsey, and M.W. Woolrich. Network mod-elling methods for fMRI. NeuroImage, 54(2):875–891, 2011.[59] M. Fiecas, H. Ombao, D. van Lunen, R. Baumgartner, A. Coimbra,and D. Feng. Quantifying temporal correlations: A test-retest eval-uation of functional connectivity in resting-state fMRI. NeuroImage,65:231–241, 2013.[60] BrainVoyager. Functional analysis: Preparation -preprocessing. Accessed: 2015-09-28.[61] M.A. Koch, D.G. Norris, and M. Hund-Georgiadis. An investigationof functional and anatomical connectivity using magnetic resonanceimaging. NeuroImage, 16(1):241–250, 2002.86Bibliography[62] M. van den Heuvel, R. Mandl, J. Luigjes, and H. Hulshoff Pol. Mi-crostructural organization of the cingulum tract and the level of de-fault mode functional connectivity. The Journal of Neuroscience,28(43):10844–10851, 2008.[63] C.J. Honey, O. Sporns, L. Cammoun, X. Gigandet, J.P. Thiran,R. Meuli, and P. Hagmann. Predicting human resting-state functionalconnectivity from structural connectivity. Proceedings of the NationalAcademy of Sciences, 106(6):2035–2040, 2009.[64] P. Skudlarski, K. Jagannathan, V.D. Calhoun, M. Hampson, B.A.Skudlarska, and G. Pearlson. Measuring brain connectivity: Diffusiontensor imaging validates resting state temporal correlations. NeuroIm-age, 43(3):554–561, 2008.[65] B. Ng, G. Varoquaux, J.B. Poline, and B. Thirion. Implications ofinconsistencies between fMRI and dMRI on multimodal connectivityestimation. In Medical Image Computing and Computer-Assisted In-tervention (MICCAI), pages 652–659. 2013.[66] K.J. Friston, A.P. Holmes, K.J. Worsley, J.-P. Poline, C.D. Frith, andR.S.J. Frackowiak. Statistical parametric maps in functional imaging:A general linear approach. Human Brain Mapping, 2(4):189–210, 1994.[67] T. Nichols and S. Hayasaka. Controlling the familywise error rate infunctional neuroimaging: a comparative review. Statistical Methodsin Medical Research, 12(5):419–446, 2003.[68] G.K. Aguirre, E. Zarahn, and M. D’Esposito. Empirical analyses ofBOLD fMRI statistics. NeuroImage, 5(3):199–212, 1997.[69] W.D. Penny, N.J. Trujillo-Barreto, and K.J. Friston. Bayesian fMRItime series analysis with spatial priors. NeuroImage, 24(2):350–362,2005.[70] L.M. Harrison, W. Penny, J. Ashburner, N. Trujillo-Barreto, andK.J. Friston. Diffusion-based spatial priors for imaging. NeuroImage,38(4):677–695, 2007.[71] C. Gossl, D.P. Auer, and L. Fahrmeir. Bayesian spatio-temporal infer-ence in functional magnetic resonance imaging. Biometrics, 57:554–562, 2001.87Bibliography[72] X. Descombes, F. Kruggel, and D.Y. von Cramon. Spatio-temporalfMRI analysis using Markov random fields. IEEE Transactions onMedical Imaging, 17(6):1028–1039, 1998.[73] W. Ou and P. Golland. From spatial regularization to anatomicalpriors in fMRI analysis. In Information Processing in Medical Imaging(IPMI), pages 88–100. 2005.[74] B. Ng, R. Abugharbieh, G. Hamarneh, and M. J. McKeown. Randomwalker based estimation and spatial analysis of probabilistic fMRI acti-vation maps. In MICCAI fMRI Data Analysis Workshop, pages 37–44.2009.[75] M.D. Fox, A.Z. Snyder, J.L. Vincent, and M.E. Raichle. Intrinsicfluctuations within cortical systems account for intertrial variabilityin human behavior. Neuron, 56(1):171–184, 2007.[76] M.D. Fox and M.E. Raichle. Spontaneous fluctuations in brain activityobserved with functional magnetic resonance imaging. Nature ReviewsNeuroscience, 8:700–711, 2007.[77] G. Deco, V.K. Jirsa, and A.R. McIntosh. Emerging concepts for thedynamical organization of resting-state activity in the brain. NatureReviews Neuroscience, 12:43–56, 2011.[78] S.M. Smith, P.T. Fox, K.L. Miller, D.C. Glahn, P.M. Fox, C.E.Mackay, N. Filippini, K.E. Watkins, R. Toro, A.R. Laird, and C.F.Beckmann. Correspondence of the brain’s functional architecture dur-ing activation and rest. Proceedings of the National Academy of Sci-ences, 106(31):13040–13045, 2009.[79] B. Ng, R. Abugharbieh, G. Varoquaux, J.-B. Poline, and B. Thirion.Connectivity-informed fMRI activation detection. In Medical Im-age Computing and Computer-Assisted Intervention (MICCAI), pages285–292. 2011.[80] K.J. Friston. Functional and effective connectivity: A review. BrainConnectivity, 1(1):13–36, 2011.[81] K. Li, L. Guo, J. Nie, G. Li, and T. Liu. Review of methods for func-tional brain connectivity detection using fMRI. Computerized MedicalImaging and Graphics, 33(2):131–139, 2009.88Bibliography[82] C. Goutte, P. Toft, E. Rostrup, F.A. Nielsen, and L.K. Hansen. Onclustering fMRI time series. NeuroImage, 9(3):298–310, 1999.[83] D. Cordes, V. Haughton, J.D. Carew, K. Arfanakis, and K. Maravilla.Hierarchical clustering to measure connectivity in fMRI resting-statedata. Magnetic Resonance Imaging, 20(4):305–317, 2002.[84] P. Golland, Y. Golland, and R. Malach. Detection of spatial activationpatterns as unsupervised segmentation of fMRI data. In Medical Im-age Computing and Computer-Assisted Intervention (MICCAI), pages110–118, 2007.[85] L. Ma, B. Wang, X. Chen, and J. Xiong. Detecting functional con-nectivity in the resting brain: a comparison between ICA and CCA.Magnetic Resonance Imaging, 25(1):47–56, 2007.[86] M.J. McKeown, S. Makeig, G.G. Brown, T.P. Jung, S.S. Kindermann,A.J. Bell, and T.J. Sejnowski. Analysis of fMRI data by blind sepa-ration into independent spatial components. Human Brain Mapping,6:160–188, 1998.[87] V. D. Calhoun, T. Adali, G. D. Pearlson, and J. J. Pekar. A method formaking group inferences from functional MRI data using independentcomponent analysis. Human Brain Mapping, 14(3):140–151, 2001.[88] S.M. Smith, K.L. Miller, S. Moeller, J. Xu, E.J. Auerbach, M.W.Woolrich, C.F. Beckmann, M. Jenkinson, J. Andersson, M.F. Glasser,D.C. Van Essen, D.A. Feinberg, E.S. Yacoub, and K. Ugurbil.Temporally-independent functional modes of spontaneous brain ac-tivity. Proceedings of the National Academy of Sciences, 109(8):3131–3136, 2012.[89] K. Wu, Y. Taki, K. Sato, Y. Sassa, K. Inoue, R. Goto, K. Okada,R. Kawashima, Y. He, A.C. Evans, and H. Fukuda. The overlappingcommunity structure of structural brain network in young healthyindividuals. PLoS ONE, 6(5):e19608, 2011.[90] V.D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre. Fastunfolding of communities in large networks. Journal of StatisticalMechanics: Theory and Experiment, 2008(10):P10008, 2008.[91] X. Yan, S. Kelley, M. Goldberg, and B.B. Biswal. Detecting overlappedfunctional clusters in resting state fMRI with connected iterative scan:89BibliographyA graph theory based clustering algorithm. Journal of NeuroscienceMethods, 199(1):108–118, 2011.[92] T. Madhyastha, Y. Cao, S. Sujitnapisatham, G. Presnyakov, W.A.Chaovalitwongse, and T. Grabowski. Link clustering to explore braindynamics using resting state functional MRI. Journal of Radiologyand Radiation Therapy, 1(2):1012, 2013.[93] Y.-Y. Ahn, J.P. Bagrow, and S. Lehmann. Link communities revealmultiscale complexity in networks. Nature, 466(7307):761–764, 2010.[94] J. Xie, S. Kelley, and B.K. Szymanski. Overlapping community detec-tion in networks: The state-of-the-art and comparative study. ACMComputing Surveys, 45(4):43, 2013.[95] R. Guimera, M. Sales-Pardo, and L.A.N. Amaral. Modularity fromfluctuations in random graphs and complex networks. Physical ReviewE, 70(2):025101, 2004.[96] A. Lancichinetti, F. Radicchi, J.J. Ramasco, and S. Fortunato. Find-ing statistically significant communities in networks. PLoS ONE,6(4):e18961, 2011.[97] C.J. Honey, J.-P. Thivierge, and O. Sporns. Can structure predictfunction in the human brain? NeuroImage, 52(3):766–776, 2010.[98] J.S. Damoiseaux and M.D. Greicius. Greater than the sum of its parts:a review of studies combining structural connectivity and resting-statefunctional connectivity. Brain Structure and Function, 213:525–533,2009.[99] B. Ng, G. Varoquaux, J.-B. Poline, and B. Thirion. A novel sparsegraphical approach for multimodal brain connectivity inference. InMedical Image Computing and Computer-Assisted Intervention (MIC-CAI), pages 707–714. 2012.[100] A. Venkataraman, Y. Rathi, M. Kubicki, C. Westin, and P. Golland.Joint modeling of anatomical and functional connectivity for popula-tion studies. IEEE Transactions on Medical Imaging, 31(2):164–182,2012.[101] L. Dodero, A. Gozzi, A. Liska, V. Murino, and D. Sona. Group-wisefunctional community detection through joint Laplacian diagonaliza-90Bibliographytion. In Medical Image Computing and Computer-Assisted Interven-tion (MICCAI), pages 708–715. 2014.[102] B. Scherrer, A. Gholipour, and S.K. Warfield. Super-resolutionin diffusion-weighted imaging. In Medical Image Computing andComputer-Assisted Intervention (MICCAI), pages 124–132. 2011.[103] A.L. Alexander, K.M. Hasan, M. Lazar, J.S. Tsuruda, and D.L.Parker. Analysis of partial volume effects in diffusion-tensor MRI.Magnetic Resonance in Medicine, 45(5):770–780, 2001.[104] S. Mori and J. Zhang. Principles of diffusion tensor imaging and itsapplications to basic neuroscience research. Neuron, 51(5):527–539,2006.[105] S. Peled and Y. Yeshurun. Superresolution in MRI: application to hu-man white matter fiber tract visualization by diffusion tensor imaging.Magnetic Resonance in Medicine, 45(1):29–35, 2001.[106] D.H.J. Poot, B. Jeurissen, Y. Bastiaensen, J. Veraart, W. Van Hecke,P.M. Parizel, and J. Sijbers. Super-resolution for multislice diffu-sion tensor imaging. Magnetic Resonance in Medicine, 69(1):103–113,2013.[107] S.N. Sotiropoulos, S. Jbabdi, J.L. Andersson, M.W. Woolrich,K. Ugurbil, and T.E.J. Behrens. RubiX: combining spatial resolu-tions for Bayesian inference of crossing fibers in diffusion MRI. IEEETransactions on Medical Imaging, 32(6):969–982, 2013.[108] V. Gupta, N. Ayache, and X. Pennec. Improving DTI resolution from asingle clinical acquisition: a statistical approach using spatial prior. InMedical Image Computing and Computer-Assisted Intervention (MIC-CAI), pages 477–484. 2013.[109] P. Coupe´, J.V. Manjo´n, M. Chamberland, M. Descoteaux, andB. Hiba. Collaborative patch-based super-resolution for diffusion-weighted images. NeuroImage, 83:245–261, 2013.[110] F. Calamante, J.-D. Tournier, G.D. Jackson, and A. Connelly. Track-density imaging (TDI): super-resolution white matter imaging usingwhole-brain track-density mapping. NeuroImage, 53(4):1233–1243,2010.91Bibliography[111] R. Zeyde, M. Elad, and M. Protter. On single image scale-up usingsparse representations. In Curves and Surfaces, pages 711–730. 2012.[112] R. Tibshirani. Regression shrinkage and selection via the lasso. Journalof the Royal Statistical Society. Series B (Methodological), pages 267–288, 1996.[113] M. Aharon, M. Elad, and A. Bruckstein. K-SVD: An algorithm fordesigning overcomplete dictionaries for sparse representation. IEEETransactions on Signal Processing, 54(11):4311–4322, 2006.[114] Y.C. Pati, R. Rezaiifar, and P. S. Krishnaprasad. Orthogonal match-ing pursuit: recursive function approximation with applications towavelet decomposition. In Asilomar Conference on Signals, Systemsand Computers, pages 40–44, 1993.[115] B.A. Landman, A.J. Huang, A. Gifford, D.S. Vikram, I.A.L. Lim,J.A.D. Farrell, J.A. Bogovic, J. Hua, M. Chen, S. Jarso, S.A. Smith,S. Joel, S. Mori, J.J. Pekar, P.B. Barker, J.L. Prince, and P.C.M. vanZijl. Multi-parametric neuroimaging reproducibility: a 3-T resourcestudy. NeuroImage, 54(4):2854 – 2866, 2011.[116] F. Deleus and M.M. Van Hulle. A connectivity-based method fordefining regions-of-interest in fMRI data. IEEE Transactions on ImageProcessing, 18(8):1760–1771, 2009.[117] R.A. Poldrack. Region of interest analysis for fMRI. Social Cognitiveand Affective Neuroscience, 2(1):67–70, 2007.[118] B. Ng, G. Hamarneh, and R. Abugharbieh. Detecting brain activationin fMRI using group random walker. In Medical Image Computing andComputer-Assisted Intervention (MICCAI), pages 331–338. 2010.[119] V. Michel, A. Gramfort, G. Varoquaux, E. Eger, C. Keribin, andB. Thirion. A supervised clustering approach for fMRI-based inferenceof brain states. Pattern Recognition, 45(6):2041–2049, 2012.[120] P.F. Neher, B. Stieltjes, I. Wolf, H.P. Meinzer, and K.H. Maier-Hein.Analysis of tractography biases introduced by anisotropic voxels. InAnnual Meeting of the International Society for Magnetic Resonancein Medicine (ISMRM). 2013.92Bibliography[121] J.V. Manjo´n, P. Coupe´, A. Buades, D.L. Collins, and M. Robles. Newmethods for MRI denoising based on sparseness and self-similarity.Medical Image Analysis, 16(1):18–27, 2012.[122] E. Garyfallidis, M. Brett, B. Amirbekian, A. Rokem, S. Van Der Walt,M. Descoteaux, and I. Nimmo-Smith. Dipy, a library for the analysisof diffusion MRI data. Frontiers in Neuroinformatics, 8(8), 2014.[123] D.C. Van Essen, S.M. Smith, D.M. Barch, T.E.J. Behrens, E. Ya-coub, and K. Ugurbil. The WU-Minn Human Connectome Project:an overview. NeuroImage, 80:62–79, 2013.[124] M.F. Glasser, S.N. Sotiropoulos, J.A. Wilson, T.S. Coalson, B. Fischl,J.L. Andersson, J. Xu, S. Jbabdi, M. Webster, and J.R Polimeni. Theminimal preprocessing pipelines for the Human Connectome Project.NeuroImage, 80:105–124, 2013.[125] Y. Behzadi, K. Restom, J. Liau, and T.T. Liu. A component basednoise correction method (CompCor) for BOLD and perfusion basedfMRI. NeuroImage, 37(1):90–101, 2007.[126] K.H. Fritzsche, P.F. Neher, I. Reicht, T. van Bruggen, C. Goch,M. Reisert, M. Nolden, S. Zelzer, H.-P. Meinzer, and B. Stielt-jes. MITK diffusion imaging. Methods of Information in Medicine,51(5):441, 2012.[127] K. Murphy, R.M. Birn, D.A. Handwerker, T.B. Jones, and P.A. Ban-dettini. The impact of global signal regression on resting state cor-relations: are anti-correlated networks introduced? NeuroImage,44(3):893–905, 2009.[128] B.P. Rogers, V.L. Morgan, A.T. Newton, and J.C. Gore. Assessingfunctional connectivity in the human brain by fMRI. Magnetic Reso-nance Imaging, 25(10):1347–1357, 2007.[129] B. Ng, G. Hamarneh, and R. Abugharbieh. Modeling brain activationin fMRI using group MRF. IEEE Transactions on Medical Imaging,31(5):1113–1123, 2012.[130] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate:a practical and powerful approach to multiple testing. Journal of theRoyal Statistical Society. Series B (Methodological), pages 289–300,1995.93Bibliography[131] Z. Liu, N. Zhang, W. Chen, and B. He. Mapping the bilateral visualintegration by EEG and fMRI. NeuroImage, 46(4):989–997, 2009.[132] L. Grady. Multilabel random walker image segmentation using priormodels. In IEEE Computer Society Conference on Computer Visionand Pattern Recognition, volume 1, pages 763–770, 2005.[133] L. Grady. Random walks for image segmentation. IEEE Transac-tions on Pattern Analysis and Machine Intelligence, 28(11):1768–1783,2006.[134] N.V. Hartvig and J.L. Jensen. Spatial mixture modeling of fMRI data.Human Brain Mapping, 11(4):233–248, 2000.[135] A.P. Dempster, N.M. Laird, and D.B. Rubin. Maximum likelihoodfrom incomplete data via the EM algorithm. Journal of the RoyalStatistical Society, Series B, 39(1):1–38, 1977.[136] K.J. Friston and W. Penny. Posterior probability maps and SPMs.NeuroImage, 19:12401249, 2003.[137] C.R. Genovese, N.A. Lazar, and T. Nichols. Thresholding of statis-tical maps in functional neuroimaging using the false discovery rate.NeuroImage, 15(4):870–878, 2002.[138] N. Toussaint, J.-C. Souplet, and P. Fillard. MedINRIA: medical imagenavigation and research tool by INRIA. In Medical Image Computingand Computer-Assisted Intervention (MICCAI), volume 7, page 280,2007.[139] L. Wang, P.D. Metzak, W.G. Honer, and T.S. Woodward. Impairedefficiency of functional networks underlying episodic memory-for-context in schizophrenia. The Journal of Neuroscience, 30(39):13171–13179, 2010.[140] H. Eichenbaum. Conscious awareness, memory and the hippocampus.Nature Neuroscience, 2:775–776, 1999.[141] K. Henke, B. Weber, S. Kneifel, H.G. Wieser, and A. Buck. Humanhippocampus associates information in memory. Proceedings of theNational Academy of Sciences, 96(10):5884–5889, 1999.[142] K. Sullivan Giovanello, D.M. Schnyer, and M. Verfaellie. A critical rolefor the anterior hippocampus in relational memory: evidence from an94BibliographyfMRI study comparing associative and item recognition. Hippocampus,14(1):5–8, 2004.[143] R.W. Heinrichs and K.K. Zakzanis. Neurocognitive deficit inschizophrenia: a quantitative review of the evidence. Neuropsychology,12(3):426, 1998.[144] F. Waters, T. Woodward, P. Allen, A. Aleman, and I. Sommer. Self-recognition deficits in schizophrenia patients with auditory halluci-nations: a meta-analysis of the literature. Schizophrenia Bulletin,38(4):741–750, 2012.[145] O. Sporns. Networks of the brain. MIT Press, 2011.[146] S. Fortunato. Community detection in graphs. Physics Reports,486(3):75–174, 2010.[147] P. Schuster and K. Sigmund. Replicator dynamics. Journal of Theo-retical Biology, 100:533–538, 1983.[148] S. Karlin. Mathematical models, problems, and controversies of evo-lutionary theory. Bulletin of the American Mathematical Society,10(2):221–273, 1984.[149] N. Meinshausen and P. Bu¨hlmann. Stability selection. Journalof the Royal Statistical Society: Series B (Statistical Methodology),72(4):417–473, 2010.[150] K. Sigmund and J. Hofbauer. The theory of evolution and dynamicalsystems. Cambridge University Press, 1988.[151] A. Torsello, S. Rota Bulo`, and M. Pelillo. Beyond partitions: Allowingoverlapping groups in pairwise clustering. In IEEE 19th InternationalConference on Pattern Recognition, pages 1–4, 2008.[152] G. Lohmann and S. Bohn. Using replicator dynamics for analyzingfMRI data of the human brain. IEEE Transactions on Medical Imag-ing, 21(5):485–492, 2002.[153] B. Thirion, P. Pinel, A. Tucholka, A. Roche, P. Ciuciu, J.-F. Mangin,and J.-B. Poline. Structural analysis of fMRI data revisited: Im-proving the sensitivity and reliability of fMRI group studies. IEEETransactions on Medical Imaging, 26(9):1256–1269, 2007.95Bibliography[154] B. Efron. The jackknife, the bootstrap and other resampling plans,volume 38. SIAM, 1982.[155] M. Lazar and A.L. Alexander. Bootstrap white matter tractography(BOOT-TRAC). NeuroImage, 24(2):524–532, 2005.[156] M. van den Heuvel, R. Mandl, and H. Hulshoff Pol. Normalized cutgroup clustering of resting-state fMRI data. PLoS ONE, 3(4):e2001,2008.[157] H.W. Kuhn. The Hungarian method for the assignment problem.Naval Research Logistics Quarterly, 2(1-2):83–97, 1955.[158] M. Rubinov and O. Sporns. Complex network measures of brain con-nectivity: uses and interpretations. NeuroImage, 52(3):1059–1069,2010.[159] Q. Wang and E. Fleury. Overlapping community structure and mod-ular overlaps in complex networks. In Mining Social Networks andSecurity Informatics, pages 15–40. 2013.[160] M.P. Van Den Heuvel and H.E.H. Pol. Exploring the brain network: areview on resting-state fMRI functional connectivity. European Neu-ropsychopharmacology, 20(8):519–534, 2010.[161] C.A. Seger. The visual corticostriatal loop through the tail of thecaudate: circuitry and function. Frontiers in Systems Neuroscience,7, 2013.[162] S.J. Palmer, J. Li, Z.J. Wang, and M.J. McKeown. Joint amplitudeand connectivity compensatory mechanisms in Parkinson’s disease.Neuroscience, 166(4):1110–1118, 2010.[163] X. Wen, L. Yao, T. Fan, X. Wu, and J. Liu. The spatial pattern of basalganglia network: A resting state fMRI study. In IEEE InternationalConference on Complex Medical Engineering, pages 43–46, 2012.[164] Y. He, Z.J. Chen, and A.C. Evans. Small-world anatomical networksin the human brain revealed by cortical thickness from MRI. CerebralCortex, 17(10):2407–2419, 2007.96


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items