Connectivity-Based Parcellation of Putamen Region usingResting State FMRIbyYiming ZhangB.Sc, Shandong University, 2013A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Electrical and Computer Engineering)The University of British Columbia(Vancouver)June 2015c© Yiming Zhang, 2015AbstractFunctional magnetic resonance imaging (fMRI) has shown great potential in study-ing the underlying neural systems. Functional connectivity measured by fMRIprovides an efficient approach to study the interactions and relationships betweendifferent brain regions. However, functional connectivity studies require accuratedefinition of brain regions, which is often difficult and may not be achieved throughanatomical landmarks. In this thesis, we present a novel framework for parcella-tion of a brain region into functional subunits based on their connectivity patternswith other reference brain regions.The proposed method takes the prior neurological information into consid-eration and aims at finding spatially continuous and functionally consistent sub-regions in a given brain region. The proposed framework relies on a sparse spa-tially regularized fused lasso regression model for feature extraction. The usuallasso model is a linear regression model commonly applied in high dimensionaldata such as fMRI signals. Compared with lasso, the proposed model further con-siders the spatial order of each voxel and thus encourages spatially and functionallyadjacent voxels to share similar regression coefficients despite of the possible spa-tial noise.In order to achieve the accurate parcellation results, we propose a process byiteratively merging voxels (groups) and tuning the parameters adaptively. In addi-tion, a Graph-Cut optimization algorithm is adopted for assigning the overlappedvoxels into separate sub-regions. With spatial information incorporated, spatiallycontinuous and functionally consistent subunits can be obtained which are desiredfor subsequent brain connectivity analysis.The simulation results demonstrate that the proposed method could reliablyiiyield spatially continuous and functionally consistent subunits. When applied toreal resting state fMRI datasets, two consistent functional subunits could be ob-tained in the putamen region for all normal subjects. Comparisons between theresults of the Parkinson’s disease group and the normal group suggest that theobtained results are in accordance with our medical assumption. The extractedfunctional subunits themselves are of great interest in studying the influence of ag-ing and a certain disease, and they may provide us deeper insights and serve as abiomarker in our future Parkinson’s disease study.iiiPrefaceThis thesis is based on one accepted conference paper:• Zhang, Yiming ; Liu, Aiping ; Tan, Sun Nee ; McKeown, Martin , Wang, Z.Jane Connectivity-Based Parcellation of Putamen using Resting State fMRIData, accepted, Proceedings of the IEEE International Symposium on Biomed-ical Imaging, 2015.The majority of this research, including literature review, model design, algo-rithm implementation, numerical simulations and draft writing, was conducted bythe author of this thesis, with suggestions from Dr. Z. Jane Wang and Dr. Mar-tin J. McKeown. Dr. Aiping Liu helped a lot on the method part in Chapter 2.Dr. Martin J. McKeown and SunNee Tan helped greatly on data collection, resultinterpretation and paper revision. Dr. Z. Jane Wang helped greatly on revisionsof the journal manuscript. Ethics approval was obtained for this research. ThefMRI data collection reported in Chapter 3.2 was covered by UBC CREB (ClinicalResearch Ethics Board) number H04-70177, with the project title as “Making theConnection: Methods to Infer Functional Connectivity in Brain Studies”.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Introduction to Popular Neural Signal Recording Techniques . . . 11.2 Introduction to FMRI and Functional Connectivity Studies . . . . 21.3 Motor Control and Parkinson’s Disease . . . . . . . . . . . . . . 81.4 Practical Motivation . . . . . . . . . . . . . . . . . . . . . . . . . 111.5 Research Objective and Methodology . . . . . . . . . . . . . . . 132 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1 Spatially Regularized Fused Lasso Model . . . . . . . . . . . . . 162.1.1 High Dimensional Linear Models . . . . . . . . . . . . . 162.1.2 Fused Lasso Method . . . . . . . . . . . . . . . . . . . . 202.1.3 Spatially Regularized Fused Lasso Method . . . . . . . . 22v2.1.4 SPG Algorithm . . . . . . . . . . . . . . . . . . . . . . . 252.2 Iterative Group Merging and Adaptive Parameter Tuning . . . . . 262.3 Graph Cut Algorithm for Assigning Overlapped Voxels . . . . . . 292.3.1 Global Minimization Problem of An Energy Function . . 302.3.2 Network Flow Problem . . . . . . . . . . . . . . . . . . . 302.3.3 Graph Cut Minimization for An Energy Function . . . . . 312.3.4 Energy Minimization in The Voxel Assignment Problem . 322.4 Extension to Regions with Three or More Sub-regions . . . . . . 343 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.1 Synthetic Datasets . . . . . . . . . . . . . . . . . . . . . . . . . . 363.2 Real FMRI Study . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2.1 Data Description . . . . . . . . . . . . . . . . . . . . . . 393.2.2 Algorithm Application . . . . . . . . . . . . . . . . . . . 393.2.3 Results on Normal Subjects . . . . . . . . . . . . . . . . 423.2.4 Results on PD Subjects . . . . . . . . . . . . . . . . . . . 423.2.5 Performance Evaluation . . . . . . . . . . . . . . . . . . 464 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . 514.1 Conclusions and Discussions . . . . . . . . . . . . . . . . . . . . 514.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54viList of TablesTable 1.1 Abbreviations of brain regions in the introduction . . . . . . . 7Table 2.1 Connectivity-based brain parcellation framework . . . . . . . . 17Table 3.1 Synthetic datasets . . . . . . . . . . . . . . . . . . . . . . . . 38Table 3.2 Best error percentage results from synthetics datasets . . . . . 38Table 3.3 Average error percentage results from synthetics Datasets . . . 39viiList of FiguresFigure 1.1 Oxyhaemoglobin and deoxyhaemoglobin in blood flow whenrest and activated. This figure is from [1] with permission. . . 3Figure 1.2 The BOLD mechanism for fMRI signals. . . . . . . . . . . . 4Figure 1.3 Example of an hemodynamic response function with simulatedfMRI signals based on [23]. . . . . . . . . . . . . . . . . . . 5Figure 1.4 Illustration of the divisions in basal ganglia and the regionsrelated with basal ganglia in the human brain. This image isfrom [13] with permission. . . . . . . . . . . . . . . . . . . . 9Figure 1.5 Illustration of two control loops of motor control in basal gan-glia in normal people. . . . . . . . . . . . . . . . . . . . . . 10Figure 1.6 Illustration of two control loops in PD patients. . . . . . . . . 11Figure 1.7 Illustration of the functional subunits in the striatum. . . . . . 13Figure 2.1 A geometrical interpretation of lasso in the left panel and ridgeregression in the right panel in the two dimensional space. . . 19Figure 2.2 A geometrical interpretation of the fused lasso model in thetwo dimensional space. . . . . . . . . . . . . . . . . . . . . . 21Figure 2.3 Illustration of a simple directed weighted graph. . . . . . . . . 31Figure 2.4 Construction of a directed weighted graph from energy mini-mization problem. . . . . . . . . . . . . . . . . . . . . . . . 32Figure 2.5 Illustration of the graph example in the voxel assignment prob-lem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 3.1 Illustration of the synthetic data generation. . . . . . . . . . . 37viiiFigure 3.2 One example of the parcellation results from syn-data 2.(a)ground truth. (b) result of the proposed spatially regularizedfused lasso algorithm. (c) result of the k-means clusteringmethod. (d) result of the modularity detection method. . . . . 40Figure 3.3 Simulation results from the spatially regularized k-means methodwith different parameter p chosen. The results are shown fromaverage error rate perspective and also minimum error rate per-spective. The upper line (blue) shows the average error ratechanging with spatial regularization parameter p, while thelower line (red) shows the minimum error rate trending. . . . 41Figure 3.4 Initial parcellation results for normal subjects N003, N004,N005 according to 3 reference regions. . . . . . . . . . . . . 43Figure 3.5 Initial parcellation results for normal subjects N007, N008,N010 according to 3 reference regions. . . . . . . . . . . . . 44Figure 3.6 Initial parcellation results for normal subjects N012, N014,N015 according to 3 reference regions. . . . . . . . . . . . . 45Figure 3.7 Final parcellation results from normal subjects N003, N004and N010. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 3.8 Final parcellation results from normal subjects N005, N007,N008,N012,N014 and N015. . . . . . . . . . . . . . . . . . . 48Figure 3.9 Final parcellation results from PD subjects P001, P002, P003,P004 ,P005 and P007. . . . . . . . . . . . . . . . . . . . . . 49Figure 3.10 Final parcellation results from PD subjects P008, P009, P010,P013 ,P015. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 3.11 The iteration process in one parcellation example (i.e., parcel-lating putamen according to SMA of the normal subject N003).The red line shows the change of the number of groups at eachiteration, the green and blue lines represent the change of thespatially regularized lasso parameter γ and the normal lassoparameter λ respectively. . . . . . . . . . . . . . . . . . . . 50ixGlossaryEEG ElectroencephalographyEMG ElectromyographyECOG ElectrocorticographyMEG MagnetoencephalographyPET Positron Emission TomographyFMRI Functional Magnetic Resonance ImagingBOLD Blood-oxygen-level dependenceHRF Hemodynamic Response FunctionICA Independent Component AnalysisSPM Statistical Parametric MappingSEM Structural Equation ModelsDCM Dynamic Causal ModelingPD Parkinson’s DiseaseBG Basal GangliaAD Alzheimer’s DiseaseROI Region-of-InterestxDLS Dorsolateral StriatumDMS Dorsomedial StriatumOF OrbitofrontalCG Cingulate GyrusSMA SensorimotorLASSO Least Absolute Shrinkage and Selection OperatorOLS Ordinary Least SquareLARS Least Angle RegressionCV Cross-ValidationSPG Smoothing Proximal GradientSNR Signal to Noise RatioxiAcknowledgmentsI would like to express my great appreciation to my supervisor, Dr. Z. Jane Wang,who has supported me throughout my master study with her constant encourage-ment and deep insight in research. I wish to thank Dr. Martin J. McKeown forhis kind financial support, precious biomedical data and enlightening suggestionsand discussions. This thesis would not have been completed or written without thegreat help and effort from Dr. Wang and Dr. McKeown.I would like to thank my examination committee members for their precioustime and constructive suggestions.Many thanks go to SunNee Tan, for her great help on data collection, biomed-ical interpretation and paper revision. I would also like to thank my dear friendsand labmates for their great help both in my life and work.At last, a thank you to my dear parents for their endless support and love. Iown my deepest gratitude to my beloved girlfriend, Youming Jiang, for sharing thehappiness as well as frustrations together with me.xiiChapter 1Introduction1.1 Introduction to Popular Neural Signal RecordingTechniquesThe human brain is the most complex organ of the body, and probably the mostcomplex thing on earth. In the history, numerous scientists and medical doctorshave devoted their whole lives to studying how the brain works from different dis-ciplines. These disciplines, including chemistry, engineering, mathematics, neu-rology, genetics, collaborate with each other to form an interdisciplinary sciencecalled neuroscience, the science of the brain. Since 1700 B.C. when the first recordabout the nervous system, Edwin Simth surgical papyrus, was written, lots of greatefforts have been made in neuroscience by many talented people in the history.Despite the progress made in neuroscience, there are still numerous mysteriesabout the brain that are unsolved. With the impressive advances of modern tech-nologies, a variety of sensors for monitoring neural activity have been developed.Neural activity recording techniques have been playing an increasingly crucial rolein neuroscience researches. Two main categories of neural activity recording tech-niques exist. One category of approaches measure neural activity by recordingelectrical activity happening in the brain or even in single neurons. These includenon-invasive methods such as Electroencephalography (EEG) and Electromyog-raphy (EMG), and also more invasive ways such as Electrocorticography (ECoG)and single-unit recording which measures the internal electrical activity of a single1neuron. Signals recorded by this category of sensors often have a high temporalresolution and low spatial resolution. The other category of approaches dependon measuring the metabolic process which will have long time constants on therecorded signals. This category of methods often fall into a broader category calledneuroimaging methods since many of them can image the function of the nervoussystem. These methods include Magnetoencephalography (MEG), Positron Emis-sion Tomography (PET) and functional Magnetic Resonance Imaging (fMRI). De-spite the relatively low temporal resolution, these type of methods often have muchhigher spatial resolution which can lead to more extracted information from thebrain.fMRI is a technique for measuring brain activity by detecting associated changesin blood flow. It has excellent spatial resolution and relatively low temporal reso-lution. fMRI is well suited for mapping human brain functions and thus has beenincreasingly employed as a popular neuroimaging tool during the past decade.1.2 Introduction to FMRI and Functional ConnectivityStudiesAround 1948, oxygen metabolism and blood flow in the brain were proved to beregulated by the brain activities with the contribution from Seymour Kety and CarlSchmidt [28]. In 1990s, Seiji Ogawa and Ken Kwong invented the fMRI techniqueand showed the potential to visualize brain activity with contrast dependent on theblood oxygen levels through the magnetic resonance imaging. The fMRI extendsthe usage of MRI to provide biological function in addition to the anatomical in-formation. The main mechanism of fMRI signals can be described as follows.When a subject performs a certain task (a cognitive task or a movement-relatedtask), some involved brain areas get activated and accordingly the oxygen con-sumption increases. After a delay of about 2 second, large increase in blood flowhappens in the local area. The increased blood flow brings more oxygen carriedby haemoglobin, so the number of oxyhaemoglobin will increase while the num-ber of deoxyhaemoglobin will consequently decrease (as shown in Figure 1.1 ).Since oxyhaemoglobin and deoxyhaemoglobin have different magnetic properties,a magnetic resonance machine can detect the change and turn it into digital signal.2Figure 1.1: Oxyhaemoglobin and deoxyhaemoglobin in blood flow when restand activated. This figure is from [1] with permission.In 2001, Logothetis confirmed that the Blood-oxygen-level dependent (BOLD)contrast mechanism reflects the neural responses triggered by stimulus [34]. Fig-ure 1.2 describes the BOLD mechanism in a flow chart. The BOLD signal of atriggered neural activity is characterized by the Hemodynamic Response Function(HRF), as seen in Figure 1.3In the early ages, fMRI was commonly used as a tool to find the activated areasin human brain when the subjects performed a certain task. During the past decade,lots of researchers have shifted their attention to studying how different brain re-gions interact, connect and coordinate with each other to perform a certain cogni-tive task rather than focusing on isolated activated brain areas. The study on theinteractions between different brain regions is referred as brain connectivity study.There are three types of brain connectivity commonly used in the brain connectiv-ity study : One is anatomical connectivity that describes the physical connectionsbetween brain sites [8]; one is effective connectivity which characterizes the causalrelation between two neural units [5]; and the third type is functional connectivitywhich is defined by temporal correlations between two brain regions [6]. Amongthese three types of brain connectivity, functional connectivity is probably the mostadopted measure in terms of brain interaction modeling. The functional connectiv-ity is often inferred by using EEG, MEG, PET and fMRI signals. Compared with3Figure 1.2: The BOLD mechanism for fMRI signals.other modalities, fMRI provides a non-invasive investigation of in vivo brain stateswith a high spatial resolution, and thus fMRI has been very popular in investigatingfunctional connectivity in the human brain.Functional connectivity modeling using fMRI signals can use either task-relatedfMRI or resting-state fMRI. Task-related fMRI requires the subject to perform aspecific task, while resting-state fMRI only needs the subjects to keep awake. Task-related fMRI is believed to reveal the corresponding brain network associated witha certain cognitive or motor task, therefore it requires a specific experimental de-sign and also brings burden to studies for aged and clinical population. On the otherhand, resting-state fMRI reveals spontaneous activity of all brain regions with theabsence of an externally prompted task [10]. Thus resting-state fMRI is preferredfor functional connectivity study in some cases.Current methods for functional connectivity analysis can be divided into twomain categories: One is the data-driven category and the other is the model-basedcategory. Data-driven methods do not rely on specific assumptions of the data40 5 10 15 20−20246810Time From Activity Onset (seconds)fMRI signalThe typical shape of an HRF HRFBaselineFigure 1.3: Example of an hemodynamic response function with simulatedfMRI signals based on [23].structure. There are two main types of approaches in the data driven category, thedecomposition methods and the clustering methods. The most widely used de-composition approach is probably Independent Component Analysis (ICA) [37].ICA decomposes the original time courses into statistically independent compo-nents (IC), the corresponding IC maps measure the correlations. By threshold-ing the IC maps, one can obtain the connectivity maps. Clustering methods forfMRI include a wide range of methods, including fuzzy clustering analysis [17],self-organizing maps [40], and neural gas network [36]. The clustering analysisnowadays is mainly considered as a preprocessing stage for connectivity studies topartition the data into different clusters that share similar connectivity patterns.Model-based methods assume certain assumptions about the data structure andare applicable to both confirmatory modeling and exploratory modeling. The lin-ear model, as a model easy to understand and easy to extend, has exerted huge5impact in brain connectivity studies. Statistical Parametric Mapping (SPM) [42],Structural Equation Models (SEM) [14] are all linear models and popular in dif-ferent scenarios of brain connectivity analysis. Dynamic Causal Modeling (DCM)[25] as a bilinear model is also popular in brain connectivity analysis. Complexnetwork analysis is another model-based method which grew very fast recentlyand draw tremendous research attentions world wide. Complex network analysis isconsidered as an advancement in connectivity analysis, since previous connectivitystudies mainly focused on how to construct a connectivity map while complex net-work analysis offers the tools to understand such a connectivity map [52]. Complexnetwork analysis, which is mathematically founded on graph theory, provides dif-ferent measures to capture the characteristics of brain connectivity networks. Oneinteresting measure is called small-worldness. A small-world network is generallydefined as a more segregated and integrated network compared with a random net-work [46]. In terms of brain connectivity map, it refers to a connectivity map thatcombines functionally specialized (segregation) sub-network with a robust num-ber of links between them (integration). Recent studies on brain connectivity net-works reported various degrees of small-worldness [2]. In addition, studies onAlzheimer’s Disease (AD) showed loss of small-worldness in AD patients [47].All these studies indicate the great potential of complex network analysis in termsof understanding the brain as connected networks.Functional connectivity analysis could be conducted at either the voxel levelor the regions-of-interest (ROI) level. Voxel-based approaches usually involve alarge number of variables. They are in general computationally inefficient, andthey are usually more sensitive to spatial noise, which could make the results un-stable. However, to accurately perform ROI-based connectivity analysis, we needto carefully define the ROIs with appropriate voxels selected. Otherwise, the sig-nal fluctuations within a single ROI may be the result of the influence of multipleunderlying functional subunits.Generally speaking, there are two ways to define the ROIs. The most straight-forward way is based on their anatomical landmarks [43]. However, the precisedefinition of a ROI is difficult and the exact boundaries of many brain regions arestill in debate. In addition, recent studies have revealed that one anatomically inde-pendent brain region could have two or more functionally separated subunits, such6Table 1.1: Abbreviations of brain regions in the introductionFull Name AbbreviationBasal Ganglia BGDorsolateral Striatum DLSDorsomedial Striatum DMSOrbitofrontal OFCingulate Gyrus CGSensori-motor SMAas the Amygdala region [9].The other way is using the data-driven approaches to define brain regions ac-cording to their functional properties, and such approaches are mainly divided intotwo categories. One is based on cluster analysis [7]. Cluster analysis often first getsconnectivity features of each voxel in areas to be parcellated using general linearregression models or Pearson pairwise correlation coefficients, and then performsclustering according to functional distances between voxels using the extracted fea-tures. Various clustering methods have been investigated in the previous studies,such as fuzzy clustering [7], k-means clustering [27], self-organized mapping [9]and maximum margin clustering [53]. However, most clustering methods demandcarefully designed denosing preprocessing. Otherwise, due to spatial noise, theyare sensitive to outliers in the data and they may fail to achieve spatially contin-uous results. The other popular category of data-driven approaches are based ongraph theory, such as the normalized cut approach [48] and the modularity detec-tion method [4]. Graph theory based approaches take the regions to be separated asa connectivity graph with each voxel representing one node in the graph and thenadopt graphical methods to separate the constructed graph. Similar to the cluster-ing methods, without the incorporation of spatial information, it is often difficultto obtain the spatially continuous subregions using graph theory methods.71.3 Motor Control and Parkinson’s DiseaseBasal Ganglia and Motor ControlBasal ganglia(BG), located deep in grey matter, is widely accepted as a primarysubstrate in the brain for motor, cognitive, and emotional processing [38]. Basalganglia contains four major divisions, and they are putamen, caudate nucleus,globus pallidus and substantia nigra. Basal ganglia is responsible for organizing themuscle-driven motor movements of the body. In another word, BG is a ‘conductor’which organizes other brain circuits to smoothly orchestrate motor behaviors [3].The illustration of the divisions in basal ganglia is shown in Figure 1.4.Among the divisions in BG, the caudate and putamen sub-regions have similarfunctions and thus together are usually taken as one area called corpus striatum orsimply striatum. Striatum is the pathway of all the inputs going to the basal gangliacircuit. Most of the inputs come from the motor cortical areas.Motor control is the process of human beings and other animals to activateand coordinate their muscles and limbs in movements through their neuromuscularsystems [45]. Motor control involves the cooperation of the motor system andthe sensory system. The motor system (i.e., the movement control system) is thepart of the central nervous system that is involved with movements. The sensorysystem (i.e., smell, touch, balance, taste, hearing and vision system) is responsiblefor processing sensory information [30]. Two control systems based on differentlearning patterns exist in motor control, both involving the participation of BG.One is goal-directed control. In goal-directed control, the behavioral selections aredetermined by relative outcome values of competing actions. Goal-directed controlis slow, effortful, serial and volitional learning. The other is habitual control. Inhabitual control, behavioral selections are determined by relative salience of learntstimulus-response associations. Habitual control is fast, parallel, and automaticlearning with low effort. According to these two control loops, the caudate andputamen integration, striatum, can be divided into two sub-divisions, dorsolateralstriatum (DLS) and dorsomedial striatum (DMS) [49]. The two control loops areillustrated in Figure 1.5.8Figure 1.4: Illustration of the divisions in basal ganglia and the regions re-lated with basal ganglia in the human brain. This image is from [13]with permission.Motor Control in Parkinson’s DiseaseParkinson’s Disease (PD) is one of the common mental disorders, and also the onewhose exact cause is still covered in the mist. PD was first introduced by andalso named after the English doctor James Parkinson in ‘An Essay on the ShakingPalsy’ in 1817 [32]. From then on, the cause, diagnosis and treatment of PD haddraw great attention from physicians, neurologists and scientists.The most obvious symptoms of PD are related with movements, including resttremor, muscle rigidity and bradykinesia. These motor symptoms are believeddue to death of dopamine producing cells in the pars compacta region of substan-9Figure 1.5: Illustration of two control loops of motor control in basal gangliain normal people.tia nigra [35], which is one part of basal ganglia. Dopamine is responsible forcarrying signals between neuros in the brain. Therefore, lack of dopamine willresult in the dysfunctionality in motor control. According to the two motor controlsystems in basal ganglia we mentioned above, PD patients are believed to sufferfrom dysfunctional habitual control since the main symptoms of PD patients canbe characterized with the loss of habitual actions such as coordinated arm swing,maintenance of balance while walking. Recent studies revealed that PD patientssuffer from dysfunctional loops of stimulus-response habit control and are trappedin goal-directed control (described in Figure 1.6) [44]. This dysfunctionality isassumed to be correlated with the loss of DLS.10Figure 1.6: Illustration of two control loops in PD patients.1.4 Practical MotivationThe main practical motivation of our study is to design a sub-division definitionmethod for DLS and DMS in the striatum and then apply the developed method toreal data to verify the hypothesis that PD patients suffer from the loss of DLS aswe mentioned above.Finding the organization of the striatum in human has always been a hot topic.Anatomically, the striatum of human is composed by the caudate and putamen, thuscannot be further divided. Recently, with the advances of noninvasive neuroimag-ing methods, the ‘functional’ organization, instead of the anatomical organizationof the striatum has been explored. These methods include functional connectivityMRI [4][20][16], T1-weighted voxel-based morphometry MRI (VBM MRI) [18]11and diffusion tensor imaging (DTI) [21][33].As we described above, functional connectivity can be ideally explored by rest-ing state fMRI. It is possible to map the striatal organization by examining theconnectivity patterns between striatum and cerebral cortex. Previous work on hu-man striatum parcellation includes [4] [20] and [16], all of them are well designedand focused on some aspects of coupling cerebral cortex with striatum. However,none of them concentrates on the two specific motor control loops that we con-cern about, i.e., the parcellation of DLS and DMS. Therefore, we design a newalgorithm for connectivity-based brain region parcellation especially suitable forour specific purpose. Then we do a pilot study with the proposed method only onthe putamen region, instead of whole striatum. The reason of doing this will beexplained below.In the motor control loops involving DLS and DMS, the putamen region mainlyreceives connections from three other brain regions in the cerebral cortex, Or-bitofrontal (OF), Cingulate Gyrus (CG) and Sensorimotor (SMA). Specifically,DLS mainly connects to SMA while DMS is mainly influenced by OF and CG.However, the spatial boundary between DLS and DMS is blurry and their exactlocations are unknown. Some voxels in DMS also receive weak connections fromSMA, and some voxels in DLS may receive weak connections from CG and OF.This scenario is illustrated in Figure 1.7. In addition, as we observed in real fMRIdata, due to the spatial signal noise, head movements and other possible artifacts,the fMRI data could be spatially corrupted with some outlier voxels. As a result,many current parcellation methods may fail to efficiently deal with such corruptedvoxels to get an accurate, spatially continuous parcellation. Therefore, we plan todesign an algorithm to parcellate the putamen region not only according to func-tional connectivity features from prior knowledge but also integrating spatial infor-mation to make it spatially continuous and functionally consistent. The proposedapproach is more flexible and reliable for the data with large spatial noise. As anadditional note, our proposed framework may not be suitable for parcellating theputamen and caudate regions simultaneously due to the incorporation of the spatialinformation. This is because the putamen and caudate are anatomically separatedand thus the proposed algorithm could be biased.12OFDMS DLSCGSMAFigure 1.7: Illustration of the functional subunits in the striatum.1.5 Research Objective and MethodologyThe technical objective of this thesis is to accurately parcellate a brain region intospatially integrated and functionally independent small sub-divisions according totheir functional interactions with other brain regions. Particularly, since parcella-tion of putamen plays a crucial role in the studies of Parkinson’s disease, we seekto find an accurate functional parcellation of putamen from which we can observethe inner organizational difference between normal people and PD patients.In order to achieve the above objective, we propose a spatially regularized con-nectivity based parcellation framework in this thesis :• Spatially Regularized Fused Lasso RegressionCommon methods of functional parcellation for brain region such as clus-tering and graph-theory methods ignore spatial information in their problemformulation. As a result, they are easily affected by irregular distributed spa-13tial noise and thus cannot promise spatially integrated results. Therefore,the main algorithm and also the first step in our framework is to acquirethe connectivity features with spatial regularization. In this thesis we em-ploy the fused lasso regression model [50] with the spatial regularizationpenalty incorporated. Lasso regression is popular in functional connectivitystudies since it fit well with high dimensional data like fMRI. Fused lassoencourages sparsity of the regression coefficients as well as sparsity of theirsuccessive differences between coefficients. In the proposed approach, weintroduce the normal lasso penalty for all voxels and the fused lasso penaltyon spatially adjacent pairs of voxels.• Iterative Group Merging and Adaptive Parameter TuningWith acquired connectivity features from spatially regularized fused lasso,voxels which share the similar connectivity patterns and neighboring siteswill be merged into groups. Then the groups are taken as new voxels andconnectivity features are re-extracted. We iteratively run this process untilwe reach groups with our desired parcellation or we reach the maximumiteration time. In addition, implementation of spatially regularized fusedlasso regression requires parameter tuning of the regularization. Ordinaryapproaches of choosing the parameter are complex and computationally in-efficient and thus can not fit in our framework. Therefore, in this thesis,we propose an adaptive parameter tuning process which tunes the parameteriteratively according to the results of each iteration.• Graph Cut Algorithm for Assigning Overlapped VoxelsAt the last step, further voxels assignment is required in overlapped areas.This is because we parcellate the brain region according to a few referenceregions, which, as a result, will lead to possible different group assignmentsin some voxels. In this thesis, we employ a graph cut method for this step.Group assignment of each overlapped voxel is decided according to its con-nectivity features as well as its spatial relationships with other voxels.The organization of the remainder thesis is as follows. In Chapter 2, the pro-posed brain region parcelltion framework is described in detail. The main algo-14rithm, the spatially regularized fused lasso model, is derived and introduced. Thenthe optimization problem is described and solved through the SPG algorithm. Next,iterative group merging and adaptive parameter tuning are explained and discussed.At last, the graph cut method is also presented. Chapter 3 reports the results whentesting the proposed framework on both synthetic datasets and real fMRI datasets.In our real fMRI study, parcellation on the putamen region was performed on ninenormal subjects and eleven PD patients. In Chapter 4, the performance of theproposed framework and results on the real fMRI study are discussed in detail.Conclusions and possible future work are given in Chapter 4.15Chapter 2MethodIn this chapter, we describe the proposed framework to separate a given brain areainto functionally consistent and spatially continuous subregions. We define thegiven brain area to be separated as the task region and other regions used to es-timate the connections with the task region as reference regions. We want to findsubsets of adjacent voxels in the task region that share similar connectivity patternsto other reference regions and iteratively merge them into groups. In the proposedalgorithm, according to prior neurological knowledge, there are several brain re-gions that share similar connectivity patterns with one functional subunit in thetask region (putamen in this thesis) so that we can obtain the functional boundarybetween the subunits.The proposed framework can be summarized in Table 2.1. We will elaborateon the individual components of the proposed framework in the following sectionsand we first describe the spatially regularized fused lasso model.2.1 Spatially Regularized Fused Lasso Model2.1.1 High Dimensional Linear ModelsThe original Least Absolute Shrinkage and Selection Operator (LASSO) algorithmis introduced mainly due to the problem in high-dimensional linear modeling.16Table 2.1: Connectivity-based brain parcellation frameworkInput : X : A T ×n matrix representing task region signal.Y1,Y2, . . . ,Ym: Each Yi is a T ×1 vector representingtime course for reference region iStep 1 : According to each reference region Yi, iterate untilfinal separation being obtained.1) Perform spatially regularized fused lasso algorithmbetween voxels in task region and one reference region.2) Merge adjacent voxels with similar connectivityweight acquired from 1) into one group.3) Treat groups from 2) as new voxels, go back to 1)Step 2 : With more than one reference region, combine allparcellation results from all reference regions byperforming graph cut algorithm.Output : Two spatially continuous functional subunits from XFirstly, think of a simple linear regression model:y = Xβ + ε (2.1)where ε ∼ N(0,σ2In), β is a p×1 column vector, y is a n× 1 column vector andX is n× p matrix indicating p variables and n observations. This linear modelhas been widely applied in detecting signals from a large dictionary during pastdecades. The original optimization problem for linear model is to minimize theordinary least square (OLS) of the model fitting to the data, so the problem couldbe formed asβˆ = argminβ‖Xβ − y‖2 (2.2)where ‖Xβ−y‖2 = (Xβ−y)T (Xβ−y). For OLS estimator, we have a closed formsolution which isβˆ = (XTX)−1XT y (2.3)The OLS estimator above works well in low dimension where usually p N.However, when p approaches N, several problems can result. One is overfittingof the model, another is collinearity among independent variables which makes17(XTX)−1 noninvertible. In addition, in high dimensional linear model, the numberof variables is often large while in most case the number of observations is small,which means n p. This makes things even worse. For example, for gene ex-pression analysis with microarray data, there could be thousands of genes but onlywith dozens of samples, since collecting one sample is quite expensive. So in highdimensional case where n p, OLS estimator is ill-posed. Fortunately, two main-stream methods exist for this situation, one is the ridge regression, another is thelasso regression.In 1970, Hoerl et al. [26] added a regularization term into the OLS estimatorto form a solution called ridge regression,βˆridgeλ= argminβ{‖Xβ − y‖2 +λ‖β‖2}(2.4)where ‖β‖2 = ∑pi=1β 2i .Ridge regression provides the trade-off between bias and variance. If one onlyconsiders about the goodness-of-fit, i.e., the bias, without constraint on β , the betacan explode and thus is susceptible to very high variance. In order to control thevariance, we might want to regularize the coefficients. Ridge regression choosesto control the L2 norm of the coefficients, in which case equation (2.4) can have aclosed form solution as follows:βˆridgeλ= (XTX +λ Ip)−1XT y (2.5)where Ip is the p-dimensional identical matrix. From the closed form solution wecan see the inclusion of the extra regularization makes the problem non-singulareven if XTX is not invertible. This is actually the original motivation for ridgeregression. Ridge regression is called as one kind of shrinkage regression techniquein which sense the shrinkage parameter λ shrinks the original coefficient size formOLS to certain extent.In 1996, Tibshirani introduced the lasso model. Similar to ridge regression,lasso also includes a regularization term for coefficients, but instead of L2 norm,lasso adopted the more aggressive L1 norm. The optimization problem for lasso is18formed as :βˆ = argminβ{‖Xβ − y‖2 +λ‖β‖1}(2.6)where ‖β‖1 =∑pi=1 |βi|. Different from ridge regression, for lasso, we often expectsparse solutions where many of the β j should be 0 in which sense lasso also doesmodel selection for us. The geometrical interpretation of difference of solutionsbetween lasso and ridge regression can be described in Figure 2.1.Figure 2.1: A geometrical interpretation of lasso in the left panel and ridgeregression in the right panel in the two dimensional space.Figure 2.1 illustrates the geometric image of lasso and ridge regression for thetwo dimensional case. x axis and y axis represent β1 and β2, the ellipse contourscorrespond to the OLS optimization function ‖Xβ − y‖2 < t1, the square in leftpanel illustrates the lasso regularization ‖β‖1 = |β1|+ |β2|< t2 while the circular inright panel shows the ridge regression regularization ‖β‖2 = β 21 +β 22 < t2. Solutionhappens when the contours first touch the square or circular. Due to the nature ofpyramid, the solution of lasso sometimes will occur at a corner in which case βiwill be equal to zero.However, unlike ridge regression, βˆ lassoλhas no closed form, due to the natureof L1 regularization. Tibshirani adopted quadratic programming techniques fromconvex optimization. In 2004, Efron et al. proposed least angle regression (LARS)[22], which computes the lasso solution efficiently.Both of ridge regression and lasso solutions are indexed by tuning parameterλ . As a result, the implementation of both methods require choosing λ . The19most commonly adopted framework is cross validation, which is widely used inthe discipline of machine learning. The most popular approach is called K-foldcross validation. In K-fold cross validation, the data will be separated into K setswith equal size. For each time, we choose the kth set as the test set and the left dataas the training set, build the model according to the training set and test the modelon the test set, calculate the prediction error ‖ fˆ (λ )−k (X)− y‖2. Then we run k timesand compute the cross-validation (CV) error:(CV Error)(λ ) = K−1K∑k=1‖ fˆ (λ )−k (X)− y‖λ2 (2.7)Then we choose the one λ with minimum CV error.In summary, in high dimensional linear regression where n p, both ridgeregression and lasso have solutions. However, lasso seems to be more attractivesince it would provide the sparse solution which reflects the inner nature in manyreal application of high dimensional linear regression.2.1.2 Fused Lasso MethodEven though lasso seems to be the most applied method in high dimensional linearregression, it has drawbacks. One of the drawbacks that people would argue is thatlasso ignores the order of the features which could be crucial in some applications.For example, in gene expression data, one would expect some genes to be highlycorrelated with each other and thus tend to put these correlated genes together in alist. In another word, the order is estimated from the data. As a result, in order totake the feature orders into consideration, Tibshirani et al. proposed a fused lassoframework [50] which adds an extra fused lasso regularization into the originallasso formulation. The fused lasso regularization is expressed as∑pj=2 |β j−β j−1| ≤t2 which penalize the differences of coefficients between each pair of features withadjacent order. With the help of Lagrange multipliers, the fused lasso optimizationproblem could be formed asβˆfused lassoλ ,γ = argminβ{‖Xβ − y‖2 +λ‖β‖1 + γp∑j=2‖Dβ‖1}(2.8)20Figure 2.2: A geometrical interpretation of the fused lasso model in the twodimensional space.where D is a p× p matrix thatD =1 −11 −1· · · · · ·1 −11 −1p×p(2.9)In equation (2.8), the first term is the original OLS estimator, to regularize thesum-of-square loss; the second term is the lasso regularization, to regularize thecoefficients and encourage coefficients to be sparse; the last term is fused lassoregularization, to encourage the sparsity in coefficients with near order.21Figure 2.2 illustrates the geometric image of fused lasso in two dimensionalcase. The ellipses represent the sum-of-square loss function, square shows regu-larization ‖β‖1 = |β1|+ |β2|< t2, the dark rectangle illustrates the fused lasso loss|β1−β2| < t2. Here we seek the first time the contours satisfy the square and therectangle. Therefore, as illustrated in Figure 2.2, sometimes the results will occurwhere the pair of the coefficients are similar. Yet, with lasso penalty, the solutionwill still be sparse.In summary, fused lasso can be classified into a more general category calledstructured lasso, since it assumed the features have a nature order. With fused lasso,one could expect a solution where sparsity exist in coefficients and also differencebetween near order coefficients. One problem in implementing fused lasso in highdimensional data is computational cost. This is especially true if one adopts K-foldcross-validation since in fused lasso there are two tuning parameters.In our application, we also concerns about the order of the coefficient, espe-cially the spacial order. Therefore, we modified the fused lasso method and adjustthe parameter tuning process according to the inner nature of the fMRI data struc-ture. In the following section, we derive our proposed spatially regularized fusedlasso method according to the properties of fMRI data.2.1.3 Spatially Regularized Fused Lasso MethodLet X = [x1,x2, · · · ,xn] be a (T × n)-dimensional data matrix with n denoting thenumber of voxels and T denoting the length of time points. X represents the fMRIsignals in the task region and xi, i = 1,2, · · · ,n represents fMRI time course ofvoxel i. Y is a (T × 1) vector representing the signals of one reference regionwhich is acquired by averaging time courses across all voxels contained in thatreference region. Let β be a (n×1) vector where each element in β represents theconnectivity coefficient between one voxel and the reference region. A standardlinear model with reference region Y is:Y = Xβ + ε (2.10)with error ε assumed to be zero mean and constant variance.The linear regression model tries to obtain the weight for each voxel in task22region in order to get the best linear fit for reference signal, i.e., smallest meansquare error. So the problem could be formed as,βˆ = argminβ‖Xβ −Y‖2 (2.11)For fMRI data, it is common that the number of voxels n in task region is largerthan the length of time points T . In this case, the problem of multicollinearitywill arise, and optimal solution of above problem will not be unique. In order toeliminate the effect of multicollinearity and make solution unique, we introducethe normal lasso penalty,n∑j=1|β j|= ‖β‖1 ≤ s1 (2.12)where s1 is a tuning parameter. This penalty will encourage sparsity in the coeffi-cient vector β , forcing small elements to be exact zero. In addition, it is well knownin fMRI analysis that similar connectivity patterns tend to exist in clusters of spa-tially adjacent voxels rather than small isolated groups of voxels [51]. However,with only lasso penalty, we cannot incorporate spatial information. Therefore, wefurther introduce spatially regularized fused lasso penalty,n∑j=1∑j<iC ji|β j−βi|< s2 (2.13)where s2 is a tuning parameter andC ji ={1, if voxel j and i are adjacent0, otherwise(2.14)Furthermore, the fused lasso penalty can be rewrote in a matrix form‖Fβ‖1 < s2 (2.15)23whereF = ( fi j) p×(p−1)2 ,p=1 −10 00 0· · · · · ·1 · · · −11 −1· · · · · ·0 0· · · · · · · · · · · ·1 −11 −11 −1p×(p−1)2 ,p(2.16)For p voxels, there are p×(p−1)2 pairs, each row in F represents one pair, so for onepair (i, j) of voxels, the non-zero elements exist if only if voxel i, j are spatiallyadjacent. In equation (2.16), each block in the matrix represents the pairs of voxeli and other voxels. The row is all zero if two voxels are not spatially adjacent.The spatially regularized fused lasso penalty comes from normal fused lassopenalty [50], and it is also a special case of graph guided fused lasso model [29].Spatially regularized fused lasso penalty encourages sparsity in difference betweenspatially adjacent voxels and thus forces adjacent voxels to belong to the samecluster.With method of Lagrange multipliers, we could form the optimization problemof spatially regularized fused lasso model as,argminβ{‖Xβ −Y‖2 +λ‖β‖1 + γ‖Fβ‖1}(2.17)where the second term encourages sparsity among elements in β , and the thirdterm encourages sparsity in differences between spatially successive elements.Note that this algorithm is suitable when n < T in fMRI analysis, because ofour assumption that voxels in the same functional subunit should share similar24connectivity weights. It makes lots of predictors suffer from collinearity whichwould allow the coefficients of predictors to be sparse and successive differencesbetween them to be sparse.Equation (2.17) leads to a quadratic programming problem with non-smoothconstraints. In order to solve this problem efficiently, we implement the SmoothingProximal Gradient (SPG) method proposed in [15].2.1.4 SPG AlgorithmThe SPG method is designed for dealing with a broad family of sparsity-inducingpenalties of complex structures. It introduces a smooth approximation to the struc-tured sparsity inducing penalty instead of optimizing the original objective functiondirectly. Additionally, it is applicable to both uni-task and multi-task sparse struc-tures. For spatially regularized fused lasso penalty, it could be taken as a generalspecial form of graph-guided-fused-lasso penalty, in which case the SPG algorithmis suitable.To implement the SPG algorithm, we first rewrite the spatially regularized lassopenalty as||Fβ ||1 = max||α||∞≤1αTFβ (2.18)where α ∈ Q = {α|‖α‖∞≤ 1} is a vector of auxiliary variables associated with‖Fβ‖1, then we could get a smooth approximation to ‖Fβ‖1 according to [14],fµ(β ) = maxα∈Q(αTFβ −µd(α)) (2.19)where µ is a positive smoothness parameter and d(α) is a smoothing functiondefined as 12‖α‖22. According to [39], fµ(β ) is smooth in β with a simple form ofthe gradient∇ fµ(β ) = FTα∗ (2.20)where α∗ is the optimal solution to the smooth approximation penalty, here weassume that for any µ , fµ(β ) is convex and continuously differentiable. For theoptimal solution α∗, we haveα∗ = S(Fβµ) (2.21)25where S is the projection operator defined as followsS(x) =x −1≤ x≤ 11 x > 1−1 x <−1(2.22)Then with smooth approximation of the L1 norm penalty, we could apply proximalgradient to solve this problem.Leth(β ) = g(β )+ fµ(β ) = ‖Y −Xβ‖2 +λ fµ(β ) (2.23)then the gradient is given as∇h(β ) = XT (Xβ −Y )+λFTα∗ (2.24)∇h(β ) is Lipschitz-continuous with Lipschitz constantL = λ (XTX)+Lµ= λmax(XTX)+‖F‖2µ(2.25)where λmax(XTX) is the largest eigenvalue of (XTX).The algorithm is formulated and described in Algorithm. 1.2.2 Iterative Group Merging and Adaptive ParameterTuningWith the lasso penalty and spatially regularized fused lasso penalty, the inferredconnectivity weights β should be sparse and also their differences sparse. In thesecond step, we merge adjacent voxels that share similar connectivity weights intogroups. To better illustrate our method, we lend some definitions from graph the-ory. We map task region as a connected graph, where each voxel is taken as onenode in graph, and edge between two nodes exists only when these two nodes areadjacent. In this network, adjacency is defined if the pair of nodes are horizontally,vertically, or diagonally adjacent in 3-dimensional space, i.e., Euclidean distancebetween two nodes is not greater than√3. We then merge those adjacency nodescorresponding to the same connectivity weights into the same group.26Algorithm 1 SPG for the Spatially Regularized Fused Lasso Penalty (modifiedfrom the algorithm in [15])Input:X, Y, F, β 0, Lipschitz constant L, desired accuracy ε , λInitializationset µ = ε2D where D = maxα∈Q12‖α‖22, θ0 = 1, w0 = βIteration:For t = 0,1,2, . . ., until convergence of β t :1. Compute ∇h(wt) = XT (Xwt −Y )+λFTα∗2. Solve the proximal operator:βt+1 = h(wt)+ 〈β −wt ,∇h(wt)〉+ L2‖β −wt‖223. Set θt+1 = 2t+3 .4. Set wt+1 = β t+1 + 1−θtθtθt+1(β t−1−β t)returnβˆ = β t+1After acquiring new groups, we average the signals contained in each groupand treat the averaged signals as new predictor variables. Then we recalculate con-nectivity weights of the updated predictor variables using the spatially regularizedfused lasso model. After merging voxels (groups), the edges in the updated graphexist when the minimal distance between voxels in two groups is not greater than√3. Repeat these processes until the number of groups reaches the expected num-ber of subunits (2 in our study) or the algorithm reaches the maximum iterationtime.However, due to the problem of tuning parameter in lasso model as well as pos-sible irregular noise of fMRI data, instead of tuning the parameter in each spatiallyregularized model, which is computationally inefficient, we adopt an adaptive pa-rameter tuning process. Our adaptive parameter tuning process tunes the parameteritself according to the results of each group merging. In detail, we first set initialregularization parameters which are relatively small. Then if the groups keep un-changed after one iteration, we will increase the penalty parameters to enforcethe regularization and relax the merging precision. Also, if after one iteration allgroups has been merged into one whole group, we will go back to the last iteration27and relax the regularization. Merging precision is controlled by a rounding param-eter δ which is defined as the precision of judging the equality of two connectiv-ity weights (For instance, merging precision δ = 1e4 means that the connectivityweights will be compared by 4 digit after decimal points).The detailed framework of parcellation with one reference region is describedin Algorithm 2.Algorithm 2 Spatially Regularized Fused Regression Based Algorithm for BrainRegion Parcellation.Input:Task region voxel level signal X , the mean time course Y of the referenceregion, the lasso penalty parameter λ , the fused lasso penalty parameter γ andthe merging precision parameter δ .Initialization:(1)Run the spatially regularized fused lasso regression model on all voxels intask region, get connectivity weights β for all voxels.(2)Round β by δ , merge adjacent voxels with same connectivity weights intoone group, get m groups. Before iteration, save original merging precisionparameter, δinitial← δ .Iteration:1: while m 6= 2 do2: Take groups as nodes, use mean value to represent groups signal. Run thespatial regularized fused lasso regression on all nodes;3: Round β by δ , merge adjacent nodes who have same connectivity weightsinto one group, get m groups.4: if m = 1 then5: γ ← γ/2, λ ← λ/2, δ ← δ ×10, redo 2 and 3;6: else7: if m remain unchanged then8: γ ← γ +1, λ ← λ +1,δ ← δ/10, redo 2 and 3;9: end if10: else11: δ ← δinitial, continue12: end if13: end while282.3 Graph Cut Algorithm for Assigning OverlappedVoxelsIn our study, several reference regions could be used for parcellation of the taskregion. Therefore, after obtaining parcellation according to each reference regionrespectively, there could be some voxels near boundary that have different groupassignments among all parcellation results. In order to further define the boundaryby combining all information of reference regions, we apply the graph cut algo-rithm to assign those voxels into groups in our framework.The graph cut algorithm is mostly applied in the field of computer vision, whereit could be used to solve image smoothing, segmentation, and many other computervision problems that could be formulated as energy minimization problems. It isbasically a combinatorial optimization technique that based on network flow theory(i.e., the max-flow min-cut theory [41]). Once an energy minimization problemcan be formulated in terms of a network flow problem, many energy functionscan be solved through the graph cut algorithm. Here we choose graph cut for thevoxel assignment problem for three main reasons. Firstly, graph cut requires theformulation of a network flow problem, in which case there exists a source nodeand a terminal node as well as usual nodes to be assigned with the label of either thesource or the terminal. This fits our framework well where we have two determinedgroups to be the source and terminal nodes and the overlapped voxels to be assignedwith labels. Secondly, the graph cut algorithm can promise the result to be spatiallyintegrated since we have to formulate the spatial relationships between voxels intothe graph cut problem formulation. At last, the graph cut algorithm is a classicapproach and easy to implement. In addition, its computational cost is little. Herewe simplify the problem with two separate regions, but it is worth to note that thisalgorithm is also applicable for three or more possible labels to assign.In the following, we first elaborate and introduce the general graph cut algo-rithm. We then formulate our problem into the graph cut framework and describethe solution.292.3.1 Global Minimization Problem of An Energy FunctionFirstly, from the very beginning, let us consider a general energy function that:E(X) =∑v∈Vgv(Xv)+∑(u,v)∈Ehuv(Xu,Xv) (2.26)where V is a set of sites, E is a set of neighboring pairs of sites, X assigns a labelto each site in V . The first term in equation (2.26) is data term while the secondis smoothing term. The problem is to find a X that minimize E(X). This energyfunction can appear in the field of different applications such as image smoothingor segmentation. Here in order to illustrate the energy function, we introduce anexample of binary image restoration.Let us consider a 2-D binary image restoration problem that we want to restorea binary image with only the noisy image Y is given. In this case, we would as-sume that the original image X is close to Y while keep the differences betweenneighboring pixels small. Therefore, to address this trade-off, we can construct anenergy function thatE(X) =∑v∈Vλ |Yv−Xv|+∑(u,v)∈Eγ|Xu−Xv| (2.27)where V is the set of all pixels and E is the set of all neighboring pairs of pixels.The main goal is to find a labeling Xv for each voxel v to minimize the energyequation above. The first term in equation (2.27) which corresponds to the dataterm in equation (2.26) regularizes the restored image X to be as close to Y aspossible. The second term which is in general the smoothing term penalizes thenear pixels to be similar. In general, minimize such a global function on all valuesis NP-Hard, thus difficult to solve. The traditional methods are all inefficient.However, in some cases, it can be solved through graph cut.2.3.2 Network Flow ProblemBefore we derive the main graph cut algorithm, we have to introduce the networkflow problem that connects the graph cut algorithm to the energy minimizationproblem.30Figure 2.3: Illustration of a simple directed weighted graph.Firstly, consider a directed weighted graph G = (V,E), where V is the set ofvertices while E is the set of edges between vertices. For a directed weighted graphG = (V,E), each edge in E has one head vertex u and one tail vertex v, the edge isdenoted as (u,v), the weight fo this edge is denoted as c(u,v). Figure 2.3 providesan example of a directed weighted graph.Then, take two vertices s and t in the directed weighted graph, and separatethe vertices into two groups S and T with respect to s and t. Also, for S and T ,V = S⋃T , S⋂T = /0, s ∈ S, t ∈ T . Such a separation is called a cut of G withrespect to s and t, denoted as (S,T ). Total sum of weights of edges from S to Tis called the cost of (S,T ). Our goal is to find a cut with minimum cost. In 1956,Ford et al. prove that the minimum cuts are given by the edges for maximum flows[24] which provides an efficient approach to find the minimum cut.2.3.3 Graph Cut Minimization for An Energy FunctionFor simplicity, here we only discuss the binary case, and consider the binary imagerestoration case. In binary image restoration, our objective is to assign all the pixels31Figure 2.4: Construction of a directed weighted graph from energy minimiza-tion problem.with 0 or 1 such that the energy function (2.27) is minimized. Therefore, constructa directed weighted graph with the source node s and the sink node t and all pixelsin the binary image. For the energy function (2.27), the graph can be constructedas in Figure 2.4.After defining the directed graph, we can see that the energy function is equalto the cut cost of this graph. Therefore, the problem of minimizing the energy istransformed into the problem of finding the minimum cut which could be given bythe edges for the maximum flow.2.3.4 Energy Minimization in The Voxel Assignment ProblemIn the proposed framework, we have two determined regions and several voxelsaround the boundary to be assigned, our goal is to assign those voxels into the twogroups. Firstly, we formalize the voxel assignment problem as the minimization ofthe following energy function,minfE( f ) =∑vDv( fv)+∑{v,q}∈NV ( fv, fq) (2.28)32Figure 2.5: Illustration of the graph example in the voxel assignment prob-lem.where f represents labels for all voxels, fv is the label of voxel v, ∑vDv( fv) is dataenergy which measures the cost of labeling voxel v as fv, ∑{v,q}∈NV ( fv, fq) is thesmooth energy that measures the extent to which the labeling f is not piecewisesmooth and N is the set of all pairs of adjacent voxels. In summary, first termin equation (2.28) tries to ensure that each voxel has been labeled correctly whilesecond term is to keep the parcellation continuity. Secondly, in order to solve theabove energy minimization problem, we have to formulate the problem as networkflow problem and define Dv( fv),V ( fv, fq). Here we construct the voxel assignmentproblem as network flow problem based on [12] and solve it using the αβ -swapalgorithm described in [11]. Figure 2.5 illustrates the graph example in voxel as-signment problem, where each small cube represents one voxel in overlapped partof the task region. α and β are labels of two determined sub-regions in task re-gion. Link ei, j between voxels exist if voxel i and voxel j are neighbors. tαi and tβirepresent links between voxel i and determined sub-region α and β respectively.In our application of graph cut, we define,Dv( fv) = 1−ρ(v, fv) (2.29)33V (α,β ) = T (α 6= β ) (2.30)where ρ(v, fv) is the Pearson correlation between voxel v and region fv, regionsignal is calculated as mean time course of the obtained subregion fv, and T (.) isPotts model.In the αβ -swap algorithm, initial label f must be given. Here we initialize thelabel fv for voxel v as α if the number of voxel v labeled as α (nα ) is larger thanthe number of that labeled as β (nβ), vice versa. If nα= nβ, we set the initial labelrandomly. Then we try to find a labeling fˆ that minimize equation (2.28) over alllabeling by swapping the labels for all voxels one by one. As described in Figure2.5, let Pαβbe the set of voxels to be assigned. There are two links connectingeach voxel with terminal α and β , noted as tαp and tβp where p indicates the voxelp. There is another possible link connecting adjacent voxels p and q together, notedas e{p,q}. According to [12], the weights assigned to the links are,tαp = Dp(α)+ ∑q∈Nv, q/∈PαβV (α, fq) for v ∈ Pαβ(2.31)tβp = Dp(β )+ ∑q∈Nv, q/∈PαβV (β , fq) for v ∈ Pαβ(2.32)e{p,q} =V (α,β ) for {v,q} ∈ N and v,q ∈ Pαβ (2.33)where Nv is the set of adjacent voxels of v, N is the set of all pairs of adjacentvoxel. Therefore, the energy function (2.28) can be solved through maximum-flowalgorithm on this graph [31].2.4 Extension to Regions with Three or MoreSub-regionsAlthough the algorithm described aforementioned is applied to our fMRI studywhere the task region is assumed to have two sub-regions, the method can be easilyextended to other brain regions with three or more sub-regions, by following somesmall modifications: First, for the parcellation of the task region according to onereference region, we only need to modify the stop criteria. The algorithm willstop merging groups when the number of groups reaches the expected number of34subregions. Secondly, graph cut is still applicable in this case. In the two sub-regions case, the αβ swap algorithm only applies on one pair of labels. While inthe case with three or more sub-regions, the αβ swap algorithm needs to run onall pairs of labels. For each pair, the process is exactly the same as we describedearlier.35Chapter 3Results3.1 Synthetic DatasetsIn this chapter, we first test the proposed method on the synthetic datasets by com-paring its performance with those of k-means clustering as well as the modularitydetection approach.In order to make simulation compliant with real fMRI conditions, we gener-ated three reference regions which connected with each other and one cubical taskregion in which one spatially continuous subunit have strong connectivity with allthree reference regions and the other part have weak connections with all referenceregions. The boundary was set to separate these two subunits.In Figure 3.1 we describe the generation of the synthetic data set. 1, 2, 3 in-dicate three connected reference regions. Y1, Y2, Y3 represent signals in referenceregions 1, 2, 3 respectively. XA represents signal in region A. XB represents signalin region B. Y and Z are source signals generated from normal distribution withzero mean and unit variance. In order to comply with real fMRI signal, we fur-ther perform temporal Gaussian smoother on the source signal. Here we set β = 1and α = 0.2. The 10×10×10 cubic task region contains 1000 voxels, and sam-ple length of each voxel was 300. Similarly, each reference region contained 300voxels with data length as 300.In the simulation, two different scenarios were considered. In the first case,data set was generated with same signal to noise ratio (SNR) for all voxels (syn-36Figure 3.1: Illustration of the synthetic data generation.data 1). Here SNR = σ2signal/σ2noise where σ2 was the variance of data. In thesecond simulation, we first randomly chose 75 voxels in each subunit of task regionin synthetic data set 1. Then we assigned those voxels with smaller SNR, i.e.,partially corrupted by outliers (syn-data 2). In order to ensure the boundary clearin synthetic data set 2, we chose the outlier voxels aside from the boundary. Resultsof two data sets were compared with those of the k-means clustering method andthe modularity detection method. Each procedure was repeated fifty times and theaveraged performance of each algorithm was evaluated by error percentage definedas,Error Percentage = # of all negative assigned voxels# of all voxels in task region(3.1)Results are shown and compared in terms of best error percentage in Table 3.2and average error percentage in Table 3.3.To implement the k-means algorithm according to [27], we first computed thePearson’s correlation coefficients between task region and each reference region atthe voxel level, and then we transformed the Pearson’s correlation ρ into Fisher’s37Table 3.1: Synthetic datasetsSNR Configurationsyn-data 1 0.5 normal spatial noisesyn-data 2 0.5(full)+0.1(part) data with few outliersTable 3.2: Best error percentage results from synthetics datasetssyn-data 1 syn-data 2Spatially Regularized Fused Lasso 0 1%K-means Clustering 0 10.90%Modularity detection 0 11.30%z statistics,z = 0.5× log 1+ρ1−ρ(3.2)Then each voxel in task region will have a corresponding feature vector where eachelement represent z score between current voxel and one voxel in reference region.k was set as 2 and Euclidean distance measure was chosen in the simulation. For themodularity detection method in [4], we constructed the undirected weighted graphby considering voxels as nodes and their pairwise Pearson’s correlation coefficientsas edge weights. Since all algorithms perform well on syn-data 1, in this thesis weonly demonstrated one representative result from syn-data 2 in Figure 3.2.In addition, to further compare with our spatially regularized framework, weperformed k-means clustering with spatial regularization. In this thesis, accord-ing to [19], we performed the spatially regularized k-means clustering method byextending the feature vector of each voxel in task region with three spatial coor-dinates, multiplied with a parameter p indicating the weights given to the spatialinformation. Figure 3.3 shows the performance of algorithm on syn-data 2 withdifferent values of parameter p, form which we could infer that with spatial infor-mation incorporated result could be improved. However, the choice of parameterp is crucial and larger value will cause loss of robustness.From the results, we noted that, without the incorporation of spatial informa-tion, it was difficult to deal with spatial noise which were not normally distributedand the outliers can not be eliminated. As a result, spatially continuous parcellation38Table 3.3: Average error percentage results from synthetics Datasetssyn-data 1 syn-data 2Spatially Regularized Fused Lasso 0.10% 2.9%K-means Clustering 0.11% 11.83%Modularity detection 0.09% 12.47 %could not be guaranteed. Nevertheless, the spatially regularized fused lasso modelwith the spatial regularization term can provide reliable parcellation with spatiallycontinuous and functionally consistent subunits.3.2 Real FMRI Study3.2.1 Data DescriptionNine healthy subjects and eleven Parkinson’s disease patients were recruited fromthe Pacific Parkinson’s Research Centre (PPRC)/Movement Disorders Clinic at theUniversity of British Columbia (UBC). During the experiment, all the subjectswere asked to lie on their back in the scanner and keep awake with their eye closed.Functional MRI data were collected from a Philips Achieva 3.0 T scanner(Philips, Best, Netherlands) equipped with a headcoil. A whole brain three-dimensionalT1-weighted image consisting of 170 axial slices with high resolution were takento facilitate the anatomical localization for each individual. BOLD contrast echo-planar (EPI) T2*-weighted images were acquired with the following specifications:repetition time 1985 ms, echo time 37 ms, flip angle 90◦, field of view (FOV)240.00mm, matrix size 128×128, and pixel size 1.9 mm×1.9 mm. Each func-tional run lasts 8 minutes. During each functional run, we obtained 36 axial sliceswith 3 mm thickness and 1 mm gap thickness. The FOV was set to contain thedorsal surface and the cerebellum of the brain.3.2.2 Algorithm ApplicationIn real fMRI data application, since left and right part of brain are symmetric, inthis pilot study, we only applied our algorithm on left part of the brain. We chose39Figure 3.2: One example of the parcellation results from syn-data 2.(a)ground truth. (b) result of the proposed spatially regularized fused lassoalgorithm. (c) result of the k-means clustering method. (d) result of themodularity detection method.OF, CG and SMA as reference regions as discussed before, trying to separate puta-men region into two functional subunits - DLS and DMS. Since we consider themotor control loops in DLS and DMS are independent, before performing anyanalysis, we first removed the global mean artifact in putamen region and all threereference regions through linear regression. Then, in order to exploit the partialcorrelation between each reference region and putamen region, we also removedartifacts from SMA in OF and CG and artifacts from CG and OF in SMA usingsame method. DLS is assumed to be functionally connected with SMA. Therefore,400 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.020.040.060.080.10.120.140.160.18Parameter pError Percentage Average Error PercentageMinimum Error PercentageFigure 3.3: Simulation results from the spatially regularized k-means methodwith different parameter p chosen. The results are shown from averageerror rate perspective and also minimum error rate perspective. Theupper line (blue) shows the average error rate changing with spatial reg-ularization parameter p, while the lower line (red) shows the minimumerror rate trending.we took SMA as reference region to define DLS in putamen region. Then for twoseparate functional subunits, in order to determine DLS, we calculated the correla-tion coefficients between each voxel contained in putamen and mean time courseof SMA. After transferring correlation coefficients to Fisher’s z-statistics, one-sideWelch’s t-test is applied on magnitude of z-statistics to test whether one region hassignificantly stronger (under 5% significance level) correlation with SMA and thisregion with stronger correlation is defined as DLS. Similarly, to define DMS inputamen region, OF and CG were selected as the reference regions according toprior knowledge. Graph cut algorithm is then essential in this case to combine twoparcellation results based on CG and OF.413.2.3 Results on Normal SubjectsFirstly, we applied our framework in parcellation of putamen in nine normal sub-jects. Figure 3.4, Figure 3.5 and Figure 3.6 demonstrate the initial parcellationresults for normal subject according to three reference regions. The axes repre-sent the 3-D coordinates. Each column shows the results from the same subjectaccording to 3 different reference regions, while each row represents results fromthe same reference region across different subjects. In each scatter plot, voxelsare shown as small circles with different colors representing groups, and the redcolor demonstrates significantly stronger correlations. It is worth to note that inour fMRI data, as a result of not using the common spatial template, the spatialcoordination of the putamen region changes from subject to subject.Figure 3.7 and Figure 3.8 show the final parcellation after the overlapping voxelassignment. In Figure 3.7 and Figure 3.8 , red points represent DMS region definedby combining initial results from CG and OF, green dots represent DLS regiondefined by initial results from SMA. Since DLS and DMS were defined by dif-ferent parcellation region, the two groups could be totally spatially separated oroverlapped. For reference, the first graph in Figure 3.7 (upper left) shows the fullputamen region in subject N003 without any parcellation. Note that although theoriginal accurate anatomical definition of putamen region may change subtly fromsubject to subject, their geometrical shapes remain almost the same.3.2.4 Results on PD SubjectsThen we applied our framework on PD subjects. Since we already illustrated theinitial parcellation results for normal subjects, here we only show the final parcel-lation results after overlapping voxels assignment. Figure 3.9 and Figure 3.10show parcellation results of putamen in eleven PD subjects. Similar to results fromnormal subject, red points represent DMS region while green dots represent DLSregion. In addition, here we use blue dots denote the overlapped voxels which ap-pears both in DLS and DMS. After we obtained the final results on PD subjects,we did quantitive analysis on the difference between results from normal and PD.42Figure 3.4: Initial parcellation results for normal subjects N003, N004, N005according to 3 reference regions.43Figure 3.5: Initial parcellation results for normal subjects N007, N008, N010according to 3 reference regions.44Figure 3.6: Initial parcellation results for normal subjects N012, N014, N015according to 3 reference regions.45Here we compared a ratio thatRatioDLS/DMS =number of voxels in DLSnumber of voxels in DMS(3.3).By comparing this ratio through student t-test, we found that the ratios in PDsubjects are significantly lower than that in normal subjects with p-value=0.0186 ,in which sense our assumption that PD suffer from loss of DLS have been verified.3.2.5 Performance EvaluationSince we have no access to the ground truth of parcellation of real fMRI data,we have to indirectly validate our results anatomically or functionally. Firstly,assuming there exists functional subunits in putamen region, then the anatomicalboundary should be stable across normal subjects. However, since our data donot use common spatial template for all subjects, spatial coordination of putamenregion changes subtly from subject to subject. As a result, it’s not easy to evaluateconsistency anatomically. The final parcellation results on normal subjects showthat DMS region (red) is consistently located on the right side of putamen regionwhile all DLS points (green) are on the left side of putamen region. Our resultsshow that the boundaries between two subunits are consistent across the subjects.In addition, since our proposed algorithm already incorporates prior knowledge (asshown in Figure 1.7) into the optimization process, the definition of DLS and DMSstrongly complies with the prior knowledge. The results demonstrate that we couldobtain spatially continuous and functionally consistent sub-regions. The proposedmethod could provide a promising way for brain parcellation.The adaptive changes of parameters of the proposed algorithm is shown inFigure 3.11 with red line representing the number of groups, green and blue linesrepresenting tuning parameters. Recall that in our algorithm, if the number ofgroups could not be reduced, we will increase the penalty parameters. It’s worth tonote that in the last few iterations, γ and λ increase dramatically since the algorithmis trying to combine those isolated outliers into the groups which demands largetunning parameters in the last several iterations.46Figure 3.7: Final parcellation results from normal subjects N003, N004 andN010.47Figure 3.8: Final parcellation results from normal subjects N005, N007,N008,N012,N014 and N015.48Figure 3.9: Final parcellation results from PD subjects P001, P002, P003,P004 ,P005 and P007.Figure 3.10: Final parcellation results from PD subjects P008, P009, P010,P013 ,P015.490 10 20 30 40 50 60 70050100150200250300350400iterations Number of groupsgammalambdaFigure 3.11: The iteration process in one parcellation example (i.e., parcel-lating putamen according to SMA of the normal subject N003). Thered line shows the change of the number of groups at each iteration, thegreen and blue lines represent the change of the spatially regularizedlasso parameter γ and the normal lasso parameter λ respectively.50Chapter 4Conclusions and Future Work4.1 Conclusions and DiscussionsIn this thesis, we present a novel framework for parcellating one brain region intofunctional subunits using fMRI data, with the basic idea of exploring functionalconnectivity patterns with other reference brain regions. The proposed approachfirst applies the spatially regularized fused lasso model to obtain functional connec-tivity relationships between voxels in the task region with other reference regions.With the learned connectivity patterns, we could group adjacent voxels with simi-lar or same connectivity weights. Then the groups containing similar connectivityadjacent voxels will be treated as new (super-)voxels and we repeat the processuntil the expected number of subregions are obtained, as described in Table 2.1. Atlast, the graph cut algorithm is applied to assign voxels into appropriate subregionsby combining the parcellation results from several reference regions.Several parameters are included in the proposed method. Parameter λ controlsthe sparsity of elements in the connectivity coefficient vector β . This parametercan be set as 0 when we only want to take sparsity in difference into consideration.It’s worth mentioning that a larger λ means higher tolerance of noise but also ahigher possibility of merging all voxels into one group. In the simulations of syn-thetic data, we set it as 1 while on real datasets we set it as 10. Another parameterγ controls the sparsity of difference between successive elements in the connec-tivity vector β . γ can not be set as 0. Otherwise, the algorithm cannot converge.51The choice of a higher initial γ could bring bias to the final results since it couldmerge too many voxels at the very first step. In our experiment, we set it as 10 forsynthetic datasets. For real fMRI datasets, due to the existence of different spa-tial noises in different subjects, we adjusted the parameter according to differentsubjects. Appropriate choice of this parameter is supported by the well-balancedgrouping results, without one group only containing less than 10 voxels or locat-ing in the middle of the other group. The third parameter is the merging precisionparameter δ , which controls the precision of merging two adjacent voxels withsimilar connectivity weights. A small δ would bring risks of merging all voxelsinto one group and a large δ would make the algorithm difficult to converge. Inour evaluations on synthetic data and real data, we both set the initial value of δas 104. It is worth noting that even though we allow adaptive adjustments of allparameters, we have to choose the initial values carefully. Inappropriate choice ofinitial parameters may lead to ill-balanced results.We have tested the proposed algorithm on synthetic data sets and comparedit with two other state-of-art methods. It is shown that the proposed algorithm iscomparable with the two other algorithms in general. In particular, the proposedmethod can well deal with the case when the spatial noise is not normally dis-tributed or when the data is corrupted with outliers. In the real fMRI study, weapplied the proposed framework to nine healthy subjects in order to define twofunctional subunits, DLS and DMS. The results are consistent with the prior med-ical knowledge. We also applied the proposed framework to the fMRI data fromeleven PD patients. Even though the results of the PD group seem to be less con-sistent than that of the normal group, our quantitative analysis shows that the PDpatients suffer from the loss of DLS. This observation is in accordance with ourassumptions, and this may serve as a potential biomarker of Parkinson’s disease.The proposed method is basically motivated by a real neurological problem,i.e., parcellating the putamen region into two functional subunits based on its con-nectivity with three other reference regions. The proposed algorithm is applicableto other neurological problems where there exists one task region to be parcellatedinto several functional subunits and a few reference regions that have different con-nectivity patterns with these functional subunits. The proposed algorithm could beapplied to both task-related fMRI data as well as the resting fMRI data. Further-52more, with suitable parameter adjustments, the proposed method does not requirethe sophisticated signal denoising procedure. However, it is worth mentioning that,although the proposed algorithm can better deal with outliers when compared withother two methods, outlier detection and removing is still essential for many con-nectivity analysis problems.4.2 Future WorkThe proposed framework can be applied to study aged people or neuro-disorderedpopulations. The extracted functional subunits themselves are of great interest instudying the influences of aging and a certain disease. Therefore, from practicalapplication perspective, we plan to apply the proposed algorithm on more data, andput more efforts to the connectivity analysis afterwards.Also, from the technical perspective, there are three possible improvements forour algorithm we could work on in the future. Firstly, since the proposed method isbased on functional connectivity, people may argue that the results for subsequentconnectivity analysis may be biased. In the future, we tends to reduce such apossible risk of bias by separating the fMRI time-series signals into two parts, anduse one part to obtain consistent parcellation results and then use the other partof data for further connectivity analysis. Secondly, since our algorithm uses meanvalue to represent the signal in reference signal, the results could be misleadingif the reference region itself contains more than one functional subunits. Thuswe plan to find a more robust representation of the reference region signals andfurther extend our framework according to this change. At last, we seek to avoidthe problem of choosing the initial parameters in our algorithm.53Bibliography[1] http://www.fmrib.ox.ac.uk/research/education. → pages viii, 3[2] S. Achard, R. Salvador, B. Whitcher, J. Suckling, and E. Bullmore. Aresilient, low-frequency, small-world human brain functional network withhighly connected association cortical hubs. The Journal of neuroscience, 26(1):63–72, 2006. → pages 6[3] G. E. Alexander and M. D. Crutcher. Functional architecture of basalganglia circuits: neural substrates of parallel processing. Trends inneurosciences, 13(7):266–271, 1990. → pages 8[4] K. A. Barnes, A. L. Cohen, J. D. Power, S. M. Nelson, Y. B. Dosenbach,F. M. Miezin, S. E. Petersen, and B. L. Schlaggar. Identifying basal gangliadivisions in individuals using resting-state functional connectivity mri.Frontiers in systems neuroscience, 4, 2010. → pages 7, 11, 12, 38[5] D. S. Bassett and E. Bullmore. Small-world brain networks. Theneuroscientist, 12(6):512–523, 2006. → pages 3[6] D. S. Bassett and E. T. Bullmore. Human brain networks in health anddisease. Current opinion in neurology, 22(4):340, 2009. → pages 3[7] R. Baumgartner, G. Scarth, C. Teichtmeister, R. Somorjai, and E. Moser.Fuzzy clustering of gradient-echo functional mri in the human visual cortex.part i: Reproducibility. Journal of Magnetic Resonance Imaging, 7(6):1094–1101, 1997. → pages 7[8] M. A. Beauchamp. An improved index of centrality. Behavioral Science, 10(2):161–163, 1965. → pages 3[9] M. Beckmann, H. Johansen-Berg, and M. F. Rushworth. Connectivity-basedparcellation of human cingulate cortex and its relation to functional54specialization. The Journal of neuroscience, 29(4):1175–1190, 2009. →pages 7[10] B. B. Biswal. Resting state fmri: a personal history. Neuroimage, 62(2):938–944, 2012. → pages 4[11] Y. Boykov and V. Kolmogorov. An experimental comparison ofmin-cut/max-flow algorithms for energy minimization in vision. PatternAnalysis and Machine Intelligence, IEEE Transactions on, 26(9):1124–1137, 2004. → pages 33[12] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy minimizationvia graph cuts. Pattern Analysis and Machine Intelligence, IEEETransactions on, 23(11):1222–1239, 2001. → pages 33, 34[13] J. W. Bradbury and S. L. Vehrencamp. Web topic 8.7, brains and decisionmaking.http://sites.sinauer.com/animalcommunication2e/chapter08.07.html.Principles of Animal Communication, Second Edition Companion Website.→ pages viii, 9[14] C. Bu¨chel and K. Friston. Modulation of connectivity in visual pathways byattention: cortical interactions evaluated with structural equation modellingand fmri. Cerebral cortex, 7(8):768–778, 1997. → pages 6[15] X. Chen, Q. Lin, S. Kim, J. G. Carbonell, E. P. Xing, et al. Smoothingproximal gradient method for general structured sparse regression. TheAnnals of Applied Statistics, 6(2):719–752, 2012. → pages 25, 27[16] E. Y. Choi, B. T. Yeo, and R. L. Buckner. The organization of the humanstriatum estimated by intrinsic functional connectivity. Journal ofneurophysiology, 108(8):2242–2263, 2012. → pages 11, 12[17] K.-H. Chuang, M.-J. Chiu, C.-C. Lin, and J.-H. Chen. Model-free functionalmri analysis using kohonen clustering neural network and fuzzy c-means.Medical Imaging, IEEE Transactions on, 18(12):1117–1128, 1999. → pages5[18] M. X. Cohen, M. V. Lombardo, and R. S. Blumenfeld. Covariance-basedsubdivision of the human striatum using t1-weighted mri. European Journalof Neuroscience, 27(6):1534–1546, 2008. → pages 1155[19] F. Deleus and M. M. Van Hulle. A connectivity-based method for definingregions-of-interest in fmri data. Image Processing, IEEE Transactions on,18(8):1760–1771, 2009. → pages 38[20] A. Di Martino, A. Scheres, D. S. Margulies, A. Kelly, L. Q. Uddin,Z. Shehzad, B. Biswal, J. R. Walters, F. X. Castellanos, and M. P. Milham.Functional connectivity of human striatum: a resting state fmri study.Cerebral cortex, 18(12):2735–2747, 2008. → pages 11, 12[21] B. Draganski, F. Kherif, S. Klo¨ppel, P. A. Cook, D. C. Alexander, G. J.Parker, R. Deichmann, J. Ashburner, and R. S. Frackowiak. Evidence forsegregated and integrative connectivity patterns in the human basal ganglia.The Journal of Neuroscience, 28(28):7143–7152, 2008. → pages 12[22] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, et al. Least angle regression.The Annals of statistics, 32(2):407–499, 2004. → pages 19[23] E. B. Erhardt, E. A. Allen, Y. Wei, T. Eichele, and V. D. Calhoun. Simtb, asimulation toolbox for fmri data under a model of spatiotemporalseparability. Neuroimage, 59(4):4160–4167, 2012. → pages viii, 5[24] L. R. Ford and D. R. Fulkerson. Maximal flow through a network. Canadianjournal of Mathematics, 8(3):399–404, 1956. → pages 31[25] K. J. Friston, L. Harrison, and W. Penny. Dynamic causal modelling.Neuroimage, 19(4):1273–1302, 2003. → pages 6[26] A. E. Hoerl and R. W. Kennard. Ridge regression: Biased estimation fornonorthogonal problems. Technometrics, 12(1):55–67, 1970. → pages 18[27] T. Kahnt, L. J. Chang, S. Q. Park, J. Heinzle, and J.-D. Haynes.Connectivity-based parcellation of the human orbitofrontal cortex. TheJournal of Neuroscience, 32(18):6240–6250, 2012. → pages 7, 37[28] S. S. Kety and C. F. Schmidt. The nitrous oxide method for the quantitativedetermination of cerebral blood flow in man: theory, procedure and normalvalues. Journal of Clinical Investigation, 27(4):476, 1948. → pages 2[29] S. Kim, K.-A. Sohn, and E. P. Xing. A multivariate regression approach toassociation analysis of a quantitative trait network. Bioinformatics, 25(12):i204–i212, 2009. → pages 24[30] B. Kolb and I. Q. Whishaw. Fundamentals of human neuropsychology.Macmillan, 2009. → pages 856[31] V. Kolmogorov and R. Zabin. What energy functions can be minimized viagraph cuts? Pattern Analysis and Machine Intelligence, IEEE Transactionson, 26(2):147–159, 2004. → pages 34[32] A. J. Lees. Unresolved issues relating to the shaking palsy on the celebrationof james parkinson’s 250th birthday. Movement Disorders, 22(S17):S327–S334, 2007. → pages 9[33] S. Lehe´ricy, M. Ducros, V. De Moortele, C. Francois, L. Thivard, C. Poupon,N. Swindale, K. Ugurbil, D.-S. Kim, et al. Diffusion tensor fiber trackingshows distinct corticostriatal circuits in humans. Annals of neurology, 55(4):522–529, 2004. → pages 12[34] N. K. Logothetis, J. Pauls, M. Augath, T. Trinath, and A. Oeltermann.Neurophysiological investigation of the basis of the fmri signal. Nature, 412(6843):150–157, 2001. → pages 3[35] J. Lotharius and P. Brundin. Pathogenesis of parkinson’s disease: dopamine,vesicles and α-synuclein. Nature Reviews Neuroscience, 3(12):932–942,2002. → pages 10[36] T. M. Martinetz, S. G. Berkovich, and K. J. Schulten. Neural-gas’ networkfor vector quantization and its application to time-series prediction. NeuralNetworks, IEEE Transactions on, 4(4):558–569, 1993. → pages 5[37] M. J. McKeown, S. Makeig, G. G. Brown, T.-P. Jung, S. S. Kindermann,A. J. Bell, and T. J. Sejnowski. Analysis of fmri data by blind separation intoindependent spatial components. Technical report, DTIC Document, 1997.→ pages 5[38] J. W. Mink. The basal ganglia: focused selection and inhibition ofcompeting motor programs. Progress in neurobiology, 50(4):381–425, 1996.→ pages 8[39] Y. Nesterov. Smooth minimization of non-smooth functions. Mathematicalprogramming, 103(1):127–152, 2005. → pages 25[40] S.-C. Ngan and X. Hu. Analysis of functional magnetic resonance imagingdata using self-organizing mapping with spatial connectivity. MagneticResonance in Medicine, 41(5):939–946, 1999. → pages 5[41] C. H. Papadimitriou and K. Steiglitz. Combinatorial optimization:algorithms and complexity. Courier Dover Publications, 1998. → pages 2957[42] W. D. Penny, K. J. Friston, J. T. Ashburner, S. J. Kiebel, and T. E. Nichols.Statistical parametric mapping: the analysis of functional brain images: theanalysis of functional brain images. Academic press, 2011. → pages 6[43] J. Rademacher, V. Caviness, H. Steinmetz, and A. Galaburda. Topographicalvariation of the human primary cortices: implications for neuroimaging,brain mapping, and neurobiology. Cerebral Cortex, 3(4):313–329, 1993. →pages 6[44] P. Redgrave, M. Rodriguez, Y. Smith, M. C. Rodriguez-Oroz, S. Lehericy,H. Bergman, Y. Agid, M. R. DeLong, and J. A. Obeso. Goal-directed andhabitual control in the basal ganglia: implications for parkinson’s disease.Nature Reviews Neuroscience, 11(11):760–772, 2010. → pages 10[45] D. A. Rosenbaum. Human motor control. Academic Press, 2009. → pages 8[46] M. Rubinov and O. Sporns. Complex network measures of brainconnectivity: uses and interpretations. Neuroimage, 52(3):1059–1069, 2010.→ pages 6[47] E. J. Sanz-Arigita, M. M. Schoonheim, J. S. Damoiseaux, S. A. Rombouts,E. Maris, F. Barkhof, P. Scheltens, and C. J. Stam. Loss ofsmall-worldnetworks in alzheimer’s disease: graph analysis of fmriresting-state functional connectivity. PloS one, 5(11):e13788, 2010. →pages 6[48] X. Shen, X. Papademetris, and R. T. Constable. Graph-theory basedparcellation of functional subunits in the brain from resting-state fmri data.Neuroimage, 50(3):1027–1035, 2010. → pages 7[49] C. A. Thorn, H. Atallah, M. Howe, and A. M. Graybiel. Differentialdynamics of activity changes in dorsolateral and dorsomedial striatal loopsduring learning. Neuron, 66(5):781–795, 2010. → pages 8[50] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, and K. Knight. Sparsity andsmoothness via the fused lasso. Journal of the Royal Statistical Society:Series B (Statistical Methodology), 67(1):91–108, 2005. → pages 14, 20, 24[51] G. Tononi, A. R. McIntosh, D. P. Russell, and G. M. Edelman. Functionalclustering: identifying strongly interactive brain regions in neuroimagingdata. Neuroimage, 7(2):133–149, 1998. → pages 2358[52] M. P. Van Den Heuvel and H. E. H. Pol. Exploring the brain network: areview on resting-state fmri functional connectivity. EuropeanNeuropsychopharmacology, 20(8):519–534, 2010. → pages 6[53] L. Wang, Q. Liu, H. Li, and D. Hu. Functional connectivity-basedparcellation of human medial frontal cortex via maximum margin clustering.In Intelligent Science and Intelligent Data Engineering, pages 306–312.Springer, 2013. → pages 759
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Connectivity-based parcellation of putamen region using...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Connectivity-based parcellation of putamen region using resting state fMRI Zhang, Yiming 2015
pdf
Page Metadata
Item Metadata
Title | Connectivity-based parcellation of putamen region using resting state fMRI |
Creator |
Zhang, Yiming |
Publisher | University of British Columbia |
Date Issued | 2015 |
Description | Functional magnetic resonance imaging (fMRI) has shown great potential in studying the underlying neural systems. Functional connectivity measured by fMRI provides an efficient approach to study the interactions and relationships between different brain regions. However, functional connectivity studies require accurate definition of brain regions, which is often difficult and may not be achieved through anatomical landmarks. In this thesis, we present a novel framework for parcellation of a brain region into functional subunits based on their connectivity patterns with other reference brain regions. The proposed method takes the prior neurological information into consideration and aims at finding spatially continuous and functionally consistent sub-regions in a given brain region. The proposed framework relies on a sparse spatially regularized fused lasso regression model for feature extraction. The usual lasso model is a linear regression model commonly applied in high dimensional data such as fMRI signals. Compared with lasso, the proposed model further considers the spatial order of each voxel and thus encourages spatially and functionally adjacent voxels to share similar regression coefficients despite of the possible spatial noise. In order to achieve the accurate parcellation results, we propose a process by iteratively merging voxels (groups) and tuning the parameters adaptively. In addition, a Graph-Cut optimization algorithm is adopted for assigning the overlapped voxels into separate sub-regions. With spatial information incorporated, spatially continuous and functionally consistent subunits can be obtained which are desired for subsequent brain connectivity analysis. The simulation results demonstrate that the proposed method could reliably yield spatially continuous and functionally consistent subunits. When applied to real resting state fMRI datasets, two consistent functional subunits could be obtained in the putamen region for all normal subjects. Comparisons between the results of the Parkinson’s disease group and the normal group suggest that the obtained results are in accordance with our medical assumption. The extracted functional subunits themselves are of great interest in studying the influence of aging and a certain disease, and they may provide us deeper insights and serve as a biomarker in our future Parkinson’s disease study. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2015-06-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0166329 |
URI | http://hdl.handle.net/2429/53970 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2015-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2015_september_zhang_yiming.pdf [ 1.25MB ]
- Metadata
- JSON: 24-1.0166329.json
- JSON-LD: 24-1.0166329-ld.json
- RDF/XML (Pretty): 24-1.0166329-rdf.xml
- RDF/JSON: 24-1.0166329-rdf.json
- Turtle: 24-1.0166329-turtle.txt
- N-Triples: 24-1.0166329-rdf-ntriples.txt
- Original Record: 24-1.0166329-source.json
- Full Text
- 24-1.0166329-fulltext.txt
- Citation
- 24-1.0166329.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0166329/manifest